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Abstract 

The  thin- shear- layer  equations  are  solved  for  a  two  dimensional 
parallel  jet  of  in  an  O^(^A)  freestream,  to  investigate  the  effect  of 
mixing  and  water  vapor  on  dissociation  in  Chemical  Oxygen- Iodine 
Lasers  (COIL) . 

Predicting  the  dissociation  of  I^  has  been  identified  as  the 
largest  source  of  error  in  COIL  performance  modeling.  Although  some 
doubt  still  remains  about  the  current  COIL  kinetics  package,  recent 
experimental  and  modeling  experiences  have  indicated  that  much  of  the 
problem  may  be  due  to  the  coupled  iodine  mixing  and  dissociation 
process.  The  numerical  code  was  used  to  study  the  combined  mixing  and 
dissociation  problem  on  a  simpler  (i.e.,  simpler  than  laser  nozzles) 
2-D,  parallel  injection  nozzle  geometry.  Comparison  of  a  set  of  1-D 
premixed  cases  and  2-D  jet  cases  with  varying  jet  velocity  ratios  (and 
therefore  mixing  rates)  indicates  that,  counter  to  the  accepted  belief, 
imperfect  mixing  which  results  in  initial  regions  of  high  I^ 
concentrai ion  leads  to  faster  dissociation  rates.  High  laser  gain, 
however,  does  require  efficient  mixing.  The  results  of  both  premixed 
and  jet-mixed  cases  with  and  without  water  vapor  in  the  0^  freestream 
demonstrate  the  strong  dependence  of  the  dissociation  efficiency  (i.e., 
the  penalty  paid  to  dissociate  I^)  on  both  the  device  geometry  and 
operating  conditions  which  challenges  the  wisdom  of  using  empirically 
determined  efficiency  factors  to  predict  COIL  performance. 
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I .  Introduction 


The  chemical  oxygen- iodine  laser  (COIL)  was  first  demonstrated  in 

1977  at  the  Air  Force  Weapons  Laboratory  (AFWL)  (1).  The  COIL  laser  is 

unique  because  it  is  a  chemically  pumped  laser  which  operates  on  an 

electronic  transition  of  the  iodine  atom  instead  of  a  vibrational  or 

rotational  transition  as  have  all  previous  chemical  lasers.  Singlet 

delta  oxygen  (0^(^A)),an  energetic  metastable  molecule,  is  created  by  a 

liquid  phase  reaction  of  chlorine  in  a  basic  hydrogen  peroxide  solution. 

Transition  of  the  O^(^A)  molecule  directly  to  the  ground  state  is  spin 

forbidden  giving  0^(V)  a  relatively  long  radiative  lifetime.  When 

mixed  with  0  (^A) ,  iodine  molecules  dissociate  into  iodine  atoms  which 
2 

are  then  pumped  rapidly  into  the  electronic  state.  Stimulated 

emission  on  the  dominant  hyperfine  transition  (F-3  — *  F-4)  of  the 

I(  - *  I(  magnetic  dipole  transition  produces  photons  ut 

1.31Sfim  making  COIL  the  shortest  wavelength  pure  chemical  laser  n 
existence.  In  1977,  the  first  COIL  laser  produced  4  mW  of  output  power. 
AFWL's  most  recent  COIL  progr2im,  Roto-COIL,  has  demonstrated  power  in 
excess  of  25  kW.  This  demonstrated  potential  for  scaling  to  high  powers 
has  generated  considerable  interest  in  COIL  devices. 

Although  the  data  base  for  designing  and  operating  COIL  devices  has 
grown  significantly  over  the  past  decade,  predicting  the  performance  of 
new  devices  or  even  modeling  the  results  of  current  devices  is  still 
somewhat  of  an  art.  Recent  experience  on  the  Roto-COIL  program 

illustrates  this  point.  The  measured  stable  resonator  power  output 

exceeded  the  predicted  performance  by  40%  even  though  the  measured 
excited  oxygen  production  was  as  expected.  The  power  available  in  the 


1 


oxygen  flow  can  be  estimated  using  the  simple  empirical  formula, 


P  -  n  (  Y  -  Y  -  N  J2  1  91  ^  (1.1) 

Aval  0  I  A  0  0  1  mole  0 

2  2  •'  2 

n^^  is  the  molar  flow  rate  of  oxygen  (all  species).  Y^  is  the  fraction 
of  oxygen  in  the  excited  state  that  is  delivered  to  the  mixing  region. 
Y^,  the  excited  oxygen  fraction  required  to  reach  lasing  threshold,  can 
be  estimated  from  the  spectroscopy  of  the  iodine  atom.  Y^  is 
approximately  0.15  at  room  temperature.  simply  the  ratio  of 

the  iodine  injected  to  the  total  oxygen  in  the  flow.  N  is  an  empirical 
value  which  is  interpreted  as  the  number  of  O^(^A)  molecules  required  to 
dissociate  one  iodine  molecule.  N  accounts  for  all  of  the  energy  lost 
during  the  mixing  and  dissociation  process  which  is  not  available  to 
pump  the  iodine  atoms.  The  lower  limit  on  N  is  set  at  2  by  the 

dissociation  energy  of  iodine.  Experimental  results  prior  to  the 

Roto-COIL  program  placed  the  value  of  N  between  5  and  7,  implying  that 
the  dissociation  was  less  than  perfectly  efficient.  With  oxygen 
production  near  the  expected  value  and  the  threshold  yield  fixed  for 
Roto-COIL  conditions,  Eq  (1.1)  implies  that  most  of  the  40%  increase  in 
power  must  be  due  to  more  efficient  iodine  dissociation,  i.e.,  N  closer 
to  2.  More  sophisticated  modeling  has  failed  to  shed  light  on  this 

apparent  discrepancy  between  the  power  available  and  and  the  measured 
extracted  power.  Even  more  troubling,  optical  measurements  of  the 

residual  in  the  cavity  indicate  that  the  1^  is  better  than  95% 
dissociated.  The  best  estimates  by  computer  models  using  the  current 
AFWL  approved  rate  package  predict  dissociation  fractions  closer  to 
50-70%. 

The  inability  to  accurately  model  the  dissociation  process  in 
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COIL  lasers  is  disturbing  as  we  begin  to  consider  larger  devices  at 
higher  power.  One  might  be  tempted  to  be  content  with  exceeding  the 
predicted  power  of  a  laser  by  40%.  Unfortunately,  laboratory  devices 
like  Roto-COIL  tend  to  to  be  designed  by  brute  force  methods.  Systems 
and  subsystems  (fluid  supply,  vacuum,  etc)  are  designed  with  sufficient 
margin  to  allow  for  some  optimization  of  the  system  performance .  As  the 
laser  is  scaled  up,  however,  the  cost  of  individual  subsystems  grows 
rapidly.  Therefore,  large  systems  are  designed  more  for  a  single 
operating  point  and  accurate  prediction  of  device  performance  becomes  an 
economic  necessity.  In  addition,  new  designs  can  not  depart  with 
confidence  from  the  current  experimental  data  base.  This  is  an 
important  consideration  as  designers  consider  trade-offs  in  H^O  trapping 
efficiency  against  transport  losses  (H^O  vapor  evolved  in  the  oxygen 
generator  rapidly  quenches  excited  iodine)  and  longer  extraction  lengths 
for  high  power  unstable  resonators. 

Because  of  these  concerns,  AFWL  initiated  a  conference  in  August 
1987  to  review  the  current  understanding  of  the  kinetic  processes 
involved  in  the  COIL  system  and  to  recommend  actions  for  improving 
modeling  capability  (2) .  While  actions  to  improve  confidence  in  several 
aspects  of  the  overall  COIL  kinetics  package  were  made,  the  summary  of 
the  COIL  kinetics  conference  concluded  that  the  largest  uncertainty  in 
modeling  COIL  lasers  was  due  to  uncertainties  in  the  iodine  dissociation 
mechanism.  One  major  problem  identified  was  the  lack  of  a  good  data 
base  to  test  the  current  dissociation  model  against.  The  kinetics  data 
base  was  assembled  from  experiments  at  conditions  quite  different  from 
laser  operating  conditions.  In  order  to  resolve  rates  of  individual 
reactions,  kinetics  experiments  are  conducted  in  low  pressure,  low 
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velocity  flow  tubes  and  where  the  mixing  considerations  are 
intentionally  designed  out,  although  it  is  not  at  all  clear  that  they 
actually  were  (c.f.  below).  The  iodine  dissociation  kinetics  are,  by 
their  nature,  heavily  influenced  by  the  mixing  process.  In  actual  laser 
nozzle  designs,  the  iodine  is  injected  transverse  to  the  primary  flow  to 
insure  rapid  entrainment  of  the  oxygen  stream  into  the  iodine  stream  and 
vise  versa.  These  flows  are  highly  three  dimensional,  making  it 
Impossible  to  sort  out  modeling  errors  in  the  solution  of  the  mixing 
problem  and  errors  due  to  the  chemical  kinetics. 

Acting  on  the  recommendations  from  the  conference ,  AFWL  began  an 
in-house  experimental  effort  to  acquire  iodine  dissociation  data  on  a 
simple,  well  diagnosed  mixing  geometry  against  which  the  current 
kinetics  package  could  be  tested.  The  goals  of  this  effort  are  to 
verify  the  existing  model  and  identify  its  weaknesses  and  develop  an 
experimental  data  base  against  which  future  changes  can  be  validated. 

This  thesis  outlines  the  development  of  a  reactive  flow  nvunerical 
model  which  will  eventually  be  used  to  aid  in  the  design  of  these 
experiments  and  to  model  the  experimental  results  although  no 
comparisons  to  experimental  data  will  be  made  here.  This  development  is 
presented  in  the  following  order.  Chapter  II  provides  a  brief  history 
of  the  development  of  the  COIL  kinetics  data  base  and  identifies  some 
of  the  sources  of  uncertainty  in  the  iodine  dissociation  process. 
Results  of  some  recent  modeling  efforts  are  presented  as  evidence  that 
mixing  plays  a  essential  role  in  the  dissociation  process.  A  brief 
description  of  the  planned  experiment  is  given.  Finally,  the  goals  the 
code  development  and  validation  effort  are  reviewed.  Chapter  III 
presents  the  detailed  development  of  the  governing  equations  used  to 
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model  the  two  dimensional,  reacting  jet  flow  field.  Chapter  IV 
discusses  the  explicit  integration  scheme  used  to  obtain  a  numerical 
solution  to  the  flow  equations.  Chapter  V  presents  the  results  of  test 
cases  used  to  verify  the  code  and  the  results  of  several  studies 
conducted  to  identify  the  roles  of  mixing  and  water  vapor  on  the  iodine 
dissociation.  Chapter  VI  summarizes  the  results  of  the  numerical  effort 
and  explores  their  Implications  to  COIL  system  performance  and  further 
experimental  and  numerical  investigations . 
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II .  Background 

Historical  Development  of  COIL  Kinetics 

In  1966,  Arnold,  Finlayson,  and  Ogryzlo  first  noted  the  near 
resonant  energy  transfer  between  electronic  states  of  molecular  oxygen 
and  atomic  iodine  (3). 

+  I(2P3^2^  - ^ 

Later,  extensive  investigation  of  this  system  by  Derwent  and  Thrush 
suggested  that  this  very  rapid  energy  transfer  (k^-  7.8x10 

(cm  /molecule  sec))  might  provide  an  excellent  pumping  mechanism  for  a 
high  energy  laser  based  on  the  dominant  hyperfine  transition  in  the 
iodine  atom,  (F-3)  — »  (F-4),  at  1.315  /m  (4)  (5).  Their  work 

quantified  the  reaction  in  Eq  (2.1)  and  showed  that  although  the  singlet 
delta  state  of  oxygen  is  nearly  resonant  with  an  electronic  state  of 
atomic  iodine,  it  is  relatively  immune  to  quenching.  The  primary  loss 
mechanism  for  02^^')  is  homogeneous  pooling. 

©2(^4)  +  02(^A)  - >  02(^S)  +  02(®E)  k  -  2.7xl0'”  (2.2) 

The  singlet  sigma  state  is  rapidly  quenched  in  the  presence  of  water. 

02(^r)  +  H2O  - »  02(^2)  +  H2O  k  -  6.7x10'^^  (2.3) 

Once  the  energy  is  transferred  to  the  iodine  atom,  the  primary 

1  * 

energy  loss  mechanisms  become  energy  pooling  of  02(  A)  and  I  , 

02(^A)  +  I*  - >  02(^2)  +  I  k  -  l.lxlO'”  (2.4) 

and  quenching  of  excited  iodine  by  water. 
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I*  +  H2O  - »  I  +  H2O  k  -  2.0x10'^^  (2.5) 

The  elimination  of  water  vapor  appears  to  be  an  important  factor  in  the 
performance  of  a  laser  system. 

Conveniently,  the  oxygen- iodine  system  also  provides  a  mechanism 
for  production  of  iodine  atoms.  Arnold  et.  al.  (3)  observed  that 
molecular  iodine  dissociates  rapidly  in  the  presence  of  excited  oxygen. 
Derwent  and  Thrush's  (4)  (5)  early  work  identified  ©2(^2)  as  the  source 
of  energy  for  this  dissociation  process.  Figure  1.  shows  the  relative 
differences  between  the  excitation  energies  of  the  two  molecules. 

02(^2)  +  I2  - >  02(^2)  +  I  +  I  k  -  4.0x10*^^  (2.6) 

©2(^2)  is  created  by  the  pooling  process  of  Eq.  (2.2)  above. 

The  remaining  requirement  to  demonstrate  lasing  was  a  source  of 

02(^A).  As  expected  for  a  near  resonant  energy  transfer,  the  reverse 

rate  for  the  pumping  reaction,  Eq  (2.1),  is  comparable  to  the  forward 

rate  (k  -  2.7x10*^^).  This  implies  that  large  number  densities  of 

O^(^A)  are  required  to  produce  a  population  inversion  in  the  iodine 

atom.  The  ratio  of  I(^P  )  to  I(^P  )  required  to  reach  threshold 

1/2  3/2 

(population  inversion)  is  determined  by  the  degeneracies  of  the  states 

(6). 

[l(^P  ) 

[l(^P  ) 

3/2 

where  [  ]  denotes  the  concentration  in  (#/cm  ) .  The  equilibrium 

condition  for  the  reaction  in  Eq  (2.1)  is  given  by 


-  -  ^  -  0.5  (2.7) 

]  ‘ 
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(2.8) 


,  [o,('i)]  [k^p.,,)]  -  [o^<=Z)]  [k'p,^^)' 


The  fraction  of  oxygen  in  the  excited  state  required  to  reach  threshold 
is  then 


0^("£)]  +  [o^<‘a)] 


0.5  k  /k 

_ r  f 

1  +  0.5  k  /k 

r'  f 


(2.9) 


The  ratio  of  the  forward  rate  to  the  reverse  rate  is  given  by  the 
equilibrium  constant  which  can  be  determined  from  statistical  mechanics. 


K  - 

•q 


k 

k 


r 


0.75 


'*oi' 


(2.10) 


At  room  temperature  (298K)  this  ratio  is  2.88.  Using  this  value,  in  Eq 
(2.9),  the  fraction  of  O^(^A)  in  the  total  oxygen  flow  required  to  lase 
a  COIL  device  must  be  at  least  0.15.  The  standard  method  of  producing 
02(^A)  for  all  of  the  early  kinetics  experiments  was  by  microwave 
discharge.  Unfortunately,  the  maximum  excited  oxygen  production  from  a 
microwave  discharge  is  only  about  0.10.  The  chemical  reaction  of 
chlorine  with  a  basic  hydrogen  peroxide  (BHP)  solution  was  found  to 
produce  large  singlet  delta  oxygen  fractions,  (s  40  -  60%)  (7), 


H2O2  +  2K0H  +  CI2  - »  2KC1  +  2H2O  +  02(^4)  (2.11) 


Although  other  methods  are  being  considered,  this  reaction  is  the  only 
easily  available  source  of  large  concentrations  of  O^(^A)  identified  to 
date . 
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Spectroscopy  of  the  Oo  -I2  -I 


Figure  1.  0,1,  and  I  Energy  Level  Diagram  for  COIL  System.  (Ref  2) 


Armed  with  this  basic  understanding  of  the  transport  and  Loss 
mechanisms  of  the  COIL  system,  in  1977  McDermott,  Pchelkin,  Benard,  and 
Bousek  designed  a  simple  longitudinal  flow  laser  which  produced  4  mW  of 
power  (1).  Within  a  year,  they  had  scaled  the  laser  to  100  W  in  a 
transverse  flow  device  (8). 

The  Iodine  Dissociation  Mechanism 

These  very  first  laser  demonstrations  provided  the  first  indication 
that  the  iodine  dissociation  mechanism  was  not  as  simple  or  as  well 
understood  as  originally  thought.  The  experiments  showed  that  high 
levels  of  I2  dissociation  were  obtained  even  in  the  presence  of  high 
concentrations  of  water  (8).  Water  is  a  strong  quencher  of  (see 

Eq  (2.3)  above).  Heidner  et.  al.  went  on  to  show  in  carefully 

controlled  flow  tube  experiments  that  the  concentration  of  62(^2)  could 
be  drastically  reduced  while  still  obtaining  significant  iodine 
dissociation  (9).  Avilds  et  al.  remeasured  the  rate  of  iodine 

dissociation  by  02(^2)  (Eq  (2.6))  and  found  it  to  be  much  smaller  than 
originally  thought  (10).  Too  small,  in  fact,  to  account  for  the 
observed  dissociation  rates. 

With  02(^21)  eliminated  as  the  primary  source  of  energy  for  I2 

dissociation,  several  teams  of  investigators  began  looking  for  an 

alternate  dissociation  mechanism  to  explain  the  experimental  results. 

It  is  clear  from  the  energy  level  diagram  in  Figure  1  that  excluding 

02(^2),  none  of  excited  states  present  have  sufficient  energy  to 

directly  dissociate  I2  by  a  single  collision.  In  1982,  Heidner  et.  al. 

proposed  a  two  step  process  involving  a  highly  vibrationally  excited  (v' 

>  43)  intermediate  of  ground  state  iodine  (denoted  as  l2^)  (H)-  Such 

1  ^ 

an  excited  state  can  be  achieved  by  collisions  with  either  02(  A)  or  I  . 
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O^Ca)  +  I2  - ^  +  I2*  k  -  7.0x10'^^  (2.12) 

I*  +  I2  - >  ^  ^2*  ^  ~  3.8x10'“  (2.13) 

Once  in  the  excited  state,  a  further  collision  with  an  02(^A)  molecule 
quickly  dissociates  the  intermediate  producing  iodine  atoms. 

I2*  +  02(^A)  - »  I  +  I  +  02(^2)  k  -  3.0x10'^°  (2.14) 

Note,  that  is  not  a  single  state  but  represents  a  number  of 

vibrational  states.  The  distribution  of  energy  in  the  vibrational 
levels  above  v"  >  43  has  been  investigated,  but  is  still  not  quantified 
(12)  (13). 

The  characteristics  of  this  dissociation  mechanism  are  quite 
different  from  the  direct  mechanism  of  Derwent  and  Thrush,  Eq  (2.6). 

4 

Pumping  to  the  intermediate  I2*  state  by  I  is  much  faster  than  by 
1  ^ 

02(  A).  Initially  however,  no  I  atoms  exist.  The  dissociation  process 
is  therefore  initially  controlled  by  the  relatively  slow  process  of 
direct  dissociation  by  02(^),  Eq  (2.6),  or  by  multiple  collisions  with 
02(^A),  Eq  (2.12)  and  Eq  (2.14).  Once  an  initial  population  of  iodine 

'ic 

atoms  is  created,  however,  they  will  be  rapidly  pumped  to  I  ,  Eq  (2.1), 

♦  * 

and  pumping  to  by  I  ,  Eq  (2.13),  will  quickly  dominate  the 

dissociation  process. 

1  1  *  . 

This  hand  off  of  control  from  02(  2)  and  02(  A)  to  I  has  a 
profound  impact  on  laser  design  and  greatly  complicates  the  task  of 
modeling  system  performance.  02(^A)  has  a  relative  long  lifetime  while 
the  lifetime  of  I  is  very  short.  Typical  laser  design  therefore 
involves  the  creation  of  02(^A)  in  a  liquid/gas  reactor,  transport 
through  cold  traps  to  remove  evolved  water  vapor,  and  injection  of 
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iodine  molecules  just  before  the  gain  region.  If  one  considered  only 
the  02(^2)  and  02(^A)  process,  the  most  efficient  mixing  scheme  would  be 
to  rapidly  mix  the  oxygen  and  iodine  streams  to  achieve  a  uniform 
mixture.  In  reality,  with  the  acceleration  of  the  dissociation  process 
due  to  1  ,  somewhat  slower  mixing  schemes  resulting  in  regions  of  high 
I2  concentrations  appear  to  result  in  faster  and  more  complete 
dissociation.  This  effect  can  be  illustrated  by  recent  modeling  and 
experimental  experiences. 

Modeling  Difficulties 

R  &  D  Associates  (RDA) ,  under  the  direction  of  AFWL,  conducted 
several  numerical  experiments  using  the  current  AFWL  approved  COIL 
kinetics  package  (14).  The  first  involved  the  comparison  of 
dissociation  results  from  a  one  dimensional  (1-D)  premixed  numerical 
code  with  the  results  of  a  two  dimensional  (2*D)  code.  The  1-D  code 
represented  an  ideal  case  of  uniform  mixing  of  the  primary  oxygen  stream 
and  the  secondary  iodine  stream.  The  initial  setup  of  the  2-D  code  is 
shown  in  Figure  2.  The  iodine  streams  are  represented  as  discrete 
stream  tubes  with  a  specified  velocity  profile.  Shear  stress  created  by 
the  velocity  differential  and  molecular  diffusion  cause  the  the  two 
streams  to  mix.  The  actual  device  involves  a  much  more  complex  three 
dimensional  (3-D)  mixing  of  sonic  transverse  jets  into  the  primary  flow 
just  upstream  of  a  supersonic  throat,  see  Figure  3.  The  2-D  calculation 
was  an  attempt  to  approximate  the  more  complicated  3-D  case.  The 
results  of  the  premixed  and  2-D  calculations  are  compared  in  Figure  4. 
The  2-D  model  predicts  a  higher  level  of  dissociation,  but  both  are  well 
below  the  experimental  value  of  nearly  100%  dissociation  by  the  nozzle 
exit  plane. 
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2.0 
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Figure  3.  Comparison  of  2-D  Mixing  Model  Geometry  to  3-D  Injection 
Scheme.  (Ref  14) 
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Fraction  of  1„  Dissociated 


X  cm 


Figure  4.  Comparison  of  1-D  Premixed  and  2-D  Mixing  Predictions  for 
Dissociation.  (Ref  14) 
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The  second  numerical  experiment  was  an  effort  to  analytically 
predict  the  results  of  Heidner's  iodine  dissociation  work  (9)  (11). 
Heidner  conducted  his  experiments  on  a  carefully  controlled,  low  speed, 
low  pressure  flow  tube.  He  used  a  microwave  discharge  to  produce  02(^A) 
molecules.  This  method  provided  a  very  clean  and  easily  controlled 
O^(^A)  source,  although  the  maximum  singlet  delta  fractions  obtainable 
were  significantly  lower  than  for  BHP/CI2  reactors.  Heidner  believed 
that  the  iodine  injectors  were  adjusted  to  remove  fluid  dynamic  and 
mixing  considerations  from  his  experiment  so  that  it  could  be  treated  as 
a  1-D  premixed  initial  condition.  RDA  modeled  the  results  of  several  of 
Heidner's  cases  using  their  1-D  premixed  code  and  the  AFV7L  COIL  kinetics 

iKr 

package  (14).  The  I  profiles  are  compared  in  Figure  5.  The  final 
levels  match  quite  well.  The  calculated  results,  however,  show  a  slower 
rise  time.  One  possible  explanation  is  that  in  the  experiment  the 
oxygen  and  iodine  were  not  well  mixed  resulting  in  local  regions  of  high 
I2  concentrations  and  gradients.  Unfortunately,  Heidner's  attempt  to 
design  a  premixed  experiment  did  not  allow  for  detailed  experimental 
investigation  of  the  initial  stages  of  the  dissociation  process  or 
precise  modeling  of  his  experimental  configuration. 

These  examples  point  out  the  need  for  increased  care  on  the  part  of 
both  experimentalist  and  modeler.  Much  of  the  early  experimental  work 
focused  solely  on  improving  the  kinetics  data  base.  Like  Heidner's 
work,  mixing  and  fluid  dynamic  considerations  were  supposedly  removed. 
Experience  demonstrates,  however,  that  the  kinetics  and  fluid  mechanics 
are  closely  coupled  making  study  of  the  kinetics  without  regard  for  the 
fluid  mechanics  both  difficult  and  unwise.  Moreover,  most  of  the 
kinetics  data  base  is  derived  from  experiments  on  flow  tubes  very 
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Figure  5.  Comparison  of  1-D  Premixed  Model  Prediction  to  Heidne: 
Experimental  Results.  (Refs  2  and  14) 


1 


similar  to  Heidner's.  Actual  laser  devices  represent  a  significant 
scaling  in  flow  velocity,  pressure,  and  number  density.  Conversely, 
nozzle  designs  for  optimum  laser  performance  usually  employ  very 
complicated  3-D  transverse  injection.  Even  simplified  representations 
of  these  geometries  have  proven  to  be  very  difficult  and  extremely 
expensive  to  accurately  model.  The  more  economical  1-D  and  2-D  models 
are  of  only  limited  value  because  they  can  be  made  to  match  experimental 
data  only  by  changing  measured  kinetics  rates.  As  a  result,  these 
methods  inspire  little  confidence  when  attempting  to  scale  current 
devices  to  higher  powers . 

A  combined  approach  which  utilizes  the  strengths  of  both  modeling 
and  experiment  is  needed  to  address  the  issues  concerning  poor  COIL 
performance  prediction.  Properly  used,  models,  analytical  or  numerical, 
are  very  powerful  investigative  tools .  Kodels  allow  an  investigator  to 
carefully  control  conditions  and  mechanisms  to  slowly  Increase  the 
complexity  of  the  problem  in  a  way  not  possible  in  experiment.  Any 
model,  however,  is  only  as  valid  as  the  physical  models  and  experimental 
data  base  upon  which  it  is  built.  Development  of  an  accurate  model  is 
essential  because  it  Implies  a  detailed  understanding  of  the  mechanisms 
responsible 

New  Approach 

A  new  combined  experimental  and  numerical  effort  has  been  initiated 
at  the  AFWL  to  once  again  address  the  iodine  dissociation  problem.  The 
approach  differs  from  the  previous  efforts.  The  objective,  at  least 
initially,  is  not  to  derive  detailed  kinetic  rate  data  but  to  tie 
modeling  employing  the  low  speed  flow  tube  kinetics  data  base  to  the 
laser  performance  data  base.  Instead  of  attempting  to  develop  a  more 
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sophisticated  numerical  scheme  to  improve  the  fluid  dynamic  modeling  of 
an  actual  laser  nozzle,  the  approach  is  to  develop  a  much  simpler,  well 
diagnosed,  experiment  which  can  be  modeled  using  existing  techniques. 
The  experimental  device  will  employ  a  BHP/Cl^  O^(^A)  generator  and 
delivery  system  to  provide  high  concentrations  of  O^(^A)  at  low  H^O 
partial  pressures.  Unlike  the  actual  laser  nozzles  however,  the  iodine 
will  be  Injected  from  a  two  dimensional,  subsonic  slit  injector. 
Figure  6  shows  a  comparison  between  the  two  geometries.  A  choke  plate 
downstream  of  the  slit  injector  insures  that  the  cavity  pressures  are 
similar  to  the  pressures  encountered  in  the  mixing  region  of  the  laser 
nozzle.  A  more  detailed  description  of  the  experiment's  design  and 
expected  operating  performance  as  they  relate  to  the  numerical  effort 
described  below  is  given  in  Appendix  A.  Admittedly,  the  two  dimensional 
mixing  scheme  will  be  less  efficient  but  the  simpler  configuration  will 
be  much  less  difficult  and  expensive  to  model.  This  approach  should 
allow  comparison  of  the  overall  characteristics  of  the  iodine 
dissociation  mechanism  to  a  numerical  simulation  when  the  fluid  dynamics 
are  well  understood.  The  results  of  this  comparison  will  allow  for  a 
better  evaluation  of  the  currently  accepted  kinetics  rate  package  and 
will  provide  a  baseline  against  which  to  test  future  changes. 

Numerical  Investigation 

This  thesis  describes  the  development  and  exercising  of  a 
two-dimensional,  reactive  flow,  thin- shear- layer  numerical  model  to 
support  the  separate  experimental  effort  described  in  Appendix  A.  The 
starting  point  for  the  numerical  model  was  a  2-D  boundary  layer  code 
developed  by  Lt  Col  Eric  Jumper  at  AFIT  (15)  (16)  to  study  HF/DF 
chemical  kinetics.  The  code  was  modified  to  solve  the  COIL  kinetics 
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Figure  6 


Comparison  of  the  2-D  and  3-D  I  Nozzle  Geometries 
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Transonic,  Transverse  Nozzle  Geometry.  (I  is 

Just  upstream  of  the  supersonic  throat).  b)  2-D, 
Jet  Nozzle  Geometry. 
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system  for  the  two  dimensional  parallel  jet  problem.  The  code  will 
eventually  be  needed  for  the  combined  experimental  and  numerical  effort, 
but  no  comparisons  to  experimental  data  will  be  made  here.  Instead, 
this  thesis  presents  the  results  of  several  cases  used  to  verify  the 
performance  of  the  code  and  to  investigate  the  roles  of  such  important 
factors  as  mixing  and  water  vapor  in  the  the  dissociation  process.  Some 
questions  are  raised  about  commonly  accepted  beliefs  about  the  iodine 
dissociation  in  COIL  devices.  Finally,  since  the  numerical  cases  were 
run  for  conditions  similar  to  those  expected  in  the  experiment,  the 
results  will  provide  some  insight  to  the  design  of  the  concurrent 
experiment. 
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III.  Governing  Equations 


Navier  Stokes  Equations 

The  derivation  of  the  differential  equations  which  govern  the  flow 
of  a  two  dimensional  plane  parallel  jet  begins  with  the  steady- state  set 
of  compressible,  multispecies  Navier-Stokes  equations  (17). 
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The  set  of  mass  diffusion  equations, 


Eq  (3.3), 


account  for  the 
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convection  and  diffusion  of  each  species  across  the  differential  control 
volume  surface  and,  in  the  case  of  reacting  flows,  the  creation  of  each 
species  due  to  chemical  reactions.  The  term  is  defined  as  the  mass 
rate  of  creation  of  species  j  per  unit  volume. 

The  energy  equation  contains  two  terms  due  to  the  presence  of 
multiple  species.  In  the  second  term  on  the  right  hand  side  of  Eq 
(3.4),  ^  is  the  dissipation  function.  This  term  accounts  for  the 
conversion  of  kinetic  energy  into  thermal  energy  by  the  action  of 
viscous  stresses.  The  third  term  is  the  change  in  enthalpy  within  the 
control  volume  due  to  diffusion  of  individual  species  across  the  control 
surface.  Note,  that  the  convection  terms  are  written  in  terms  of  the 
total  enthalpy  of  the  flow  where  ^  “  Zj  conversion  of 
"chemical"  energy  to  thermal  energy  is,  by  definition,  included  in  these 
terms.  The  energy  equation  can  be  rewritten  employing  the  mass 
diffusion  equations  to  express  the  addition  of  chemical  energy  to  the 
flow  in  a  more  explicit  form. 

Substituting  for  the  enthalpy  in  terms  of  the  enthalpies  of  each 
species,  and  expanding  both  the  convection  and  species  diffusion  terms 
using  the  chain  rule  gives. 
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The  last  term  in  Eq  (3.6)  can  be  solved  for  from  the  mass  diffusion 
equation  (3.3). 
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7  ^  ^ 


(3.7) 


The  last  term  in  Eq  (3.7)  is  the  change  in  enthalpy  due  to  chemical 
reactions.  Eq  (3.7)  can  be  substituted  into  Eq  (3.6)  for  further 
simplification.  Finally,  if  each  component  of  the  mixture  is  assumed  to 
behave  like  an  ideal  gas,  then,  by  definition  of  the  specific  heat, 


cp. 


and 


cp  -  Z  W  cp 
J 


The  energy  equation  becomes ; 


pcp 


3P 

1 


t::  + 


2  W  h 
.  J  J 


dV  dh  ] 

ax^  ax^  I 

1  I  •' 


+ 


(3.8) 


Thin- Shear -Laver  Equations 

The  complexity  of  the  Navier-Stokes  equations  can  be  significantly 
reduced  for  a  commonly  occurring  class  of  flows  in  which  the  viscous 
forces  are  confined  to  a  thin  layer  relative  to  the  overall  flow  field. 
This  class  is  referred  to  as  the  boundary- layer  approximation  in  the 
case  of  immersed  bodies  or,  more  generally,  as  the  thin-shear- layers 
approximation  in  the  case  of  jets,  wakes,  and  mixing-layers  (18).  The 
plane  parallel  jet  was  in  fact  selected  as  the  mixing  scheme  to  be  used 
in  this  investigation  because  it  has  been  extensively  studied  (19)  (20) , 
and  because  it  can  be  modeled  using  the  thin- shear- layer  equations. 

The  reduced  set  of  governing  equations  for  thin- shear -layers  is 
derived  by  placing  the  following  constraints  on  the  Navier-Stokes 
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equations.  First,  the  viscous -layers  (including  the  thermal- layers  and 
diffusion- layers)  are  thin  relative  to  the  characteristic  streamwise 
dimension  of  the  flow  field.  Let  S  be  thicknesses  of  the  viscous 


shear- layer  and  L  be  a  characteristic  streamwise  length,  then  this  first 
constraint  implies  that  6/L  «  1  .  Similar  requirements  are  placed  on 
the  thermal  and  diffusion  boundary- layers .  The  second  requirement  is 
that  the  largest  viscous  term  be  of  the  same  order  of  magnitude  as  the 
convective  terms.  The  thin- shear- layer  approximations  can  be  applied  to 
the  Navier-Stokes  equations  in  a  number  of  ways.  An  order  of  magnitude 
analysis  is  used  here  on  a  term  by  term  basis  (See  Schlichting  (18)  for 
more  details).  5/L  «  1  implies  that  the  streamwise  dimension,  x,  is 
0(L)  while  the  transverse  dimension,  y,  is  0(5).  The  streamwise 
velocity  component  must  be  0(1).  The  terms  of  the  continuity  equation 
must  be  of  the  same  order.  3u/5x  a  0(1)  so  the  transverse  velocity 
component,  v,  must  be  of  0(5/L).  The  pressure  gradient  term,  3P/3x,  in 
the  streamwise  momentum  equation  must  be  balance  the  convection  terms  in 
the  inviscid  limit,  this  implies  that  3P/3x  =  0(1).  The  pressure 
gradient  term  in  transverse  momentum  equation  can  not  be  of  higher  order 
than  any  of  the  other  terms  which  is  0(5/L).  The  streamwise  pressure 
gradient  can  then  be  rewritten  as  dP/dx.  In  general,  for  large  Reynolds 
numbers,  the  terms  3^(  )/3x^  are  of  0(5/L)  or  smaller.  Neglecting 
terms  of  0(5/L)  and  smaller  results  in  the  thin- shear- layer  equations. 

I-  (pu)  +  I-  (pv)  -  0  (3.9) 


pu 


3u 

3x 


3v 

pv  -r- 

3y 


dx  3y 


(3.10) 
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(3.11) 


dw 


3W 


-  +  W  + 


-  f 

dy  [ 


pD. 


aw 

ay^ 


(3.12) 


von  Mises  Transformation 

The  2-D,  steady,  thin-shear- layer  equations  can  be  recast  into  a 
more  convenient  form  for  numerical  integration  using  the  compressible 
stream  function,  ^  .  The  stream  function  is  defined  as  the  mass  flow 
(per  unit  depth)  between  the  origin  and  any  point.  By  definition,  there 
is  no  net  flow  of  mass  across  lines  of  constant  r/>  .  This  does  not, 
however,  rule  out  net  transport  of  individual  species.  The 
transformation  from  cartesian  coordinates,  x-y,  to  stream  function 
coordinates,  x-^,  is  accomplished  using  che  von  Mises  transformation 
(15)  (16). 


drf) 


Then  for  any  arbitrary  function 


af  _  af  ^ 
ax  a^  ax 


and 


ay 


pu 


f  -  f(x,y), 

af  ax  _  _  3f  ^ 
ax  ax  a^  ax 


af_afa^  afax_  ^ 
ay  a0  ay  ^  ax  ay  di> 

Applying  the  von  Mises  transformation  to  the  boundary -layer  equations 
results  in, 
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(3.13) 


dx  pu  dx  ^  di> 


aw. 

ax-^ 


pu  j 


( 


pD  pu 
J 


aw  \ 

di/>'  j 


(3.14) 


3T 

ax 


1  ^  ^  /££U 

pcp  dx  Cp 


f  ^ 

[  dr/> 


j 


Thermodynamic  Properties 

The  density  of  the  mixture  of  gases  was  calculated  from  the  ideal 
gas  equation, 


P 


P 

R  T  W 

mix 


where  H  is  the  molecular  weight  of  the  mixture, 
mix 


(3.16) 


M  -  y  y  W 

mix  ^  ''i  i 
i 

The  specific  heat  and  enthalpy  for  each  of  the  component  species  of 
the  gas  were  calculated  using  fourth  order  polynomial  fits  to  the  data 
tabulated  in  the  JANAF  Thermodynamic  Tables  (21). 

Cp  -  a  +  a  T  +  a  T^  +  a  T^  +  a  T*  (3.17) 

i  0  1  2  3  4 

h  -  b  +  bT  +  bT^  +  bT^  +  bT*  +  E  (3.18) 

1012311 

The  units  from  the  JANAF  tables  are  cp^  [-]  (kcal/mole  K)  and 
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[“]  (kcal/mole) .  The  enthalpy  data  is  referenced  at  298.15  K.  is 

the  heat  of  fonnation  of  species  i  (22), 

E  -  AH  +  hcN  1/ 

1  298  A  i 


where 

AH°  -  heat  of  formation  of  the  ground  state 
f298  ° 

-  excitation  energy  of  the  i^^  excited  state  (cm 

The  excitation  energies  for  the  O^(^A)  and  O^(^S)  states  of  oxygen  are 

7,882  cm  ^  and  13,121  cm  respectively  (20).  The  excitation  energy  of 

*  -1 

the  excited  iodine  atom,  I  ,  is  7,603  cm  As  stated  in  Chapter  II, 

however,  the  intermediate  excited  state  of  the  iodine  molecule,  I  ^ , 

2 

actually  represents  a  number  of  vibrational  states  with  v"  >  43.  Little 
is  known  about  the  distribution  of  energies  for  the  I^^  states. 
Inspection  of  Eqs  (2.8)  and  (2.9)  indicates  that  the  excitation  energy 
should  be  approximately  the  same  as  that  of  O^(^A)  or  I*.  The 

excitation  energy  for  I^*  is  therefore  set  equal  to  that  for  I  ,  7,603 

-1 

cm 

The  specific  heat  of  the  mixture  of  gases  is 

cp  -  Z 

1 

The  enthalpy  of  species  i  per  gram  is,  h  /M^.  The  specific  heat  and 
enthalpy  are  converted  to  the  CCS  units  (cm^/s^K)  and  (cmVs^) 
respectively. 

Transport  Coefficients 

The  transport  coefficients  for  viscous,  thermal,  and  mass  diffusion 
are  derived  form  the  kinetic  theory  of  gases  (24).  Although  expressions 
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for  the  transport  coefficients  can  be  obtained  by  modeling  the  molecules 
of  a  gas  as  rigid  spheres,  the  correlation  of  the  results  with 
experimental  data  is  poor.  Better  correlation  is  obtained  using  the 
Chapman- Enskog  theory  which  derives  the  transport  coefficients  in  terms 
of  the  potential  energy  of  interaction  between  two  molecules.  The 
potential  energy,  <p,  is  related  to  the  force  of  interaction  between 
molecules  by  F  -  -  dep/dr  where  r  is  the  intermolecular  distance. 
Kinetic  theory  does  not  provide  a  formal  expression  for  ip(r)  however, 
there  is  an  extensive  empirical  data  base  in  terms  of  the  Lennard- Jones 
potential . 

<7  is  a  characteristic  diameter  and  e  is  the  maximvim  energy  of 
attraction.  Note  that  this  empirical  function  displays  the  expected 
characteristics  of  weak  attraction  at  large  separations  and  and  strong 
repulsion  at  short  separations. 

Viscosity.  The  expression  for  the  viscosity  of  component  i  from 
the  Chapman- Enskog  theory  is 


i 


2.6693  X  10 


-5 


n 


(3.19) 


The  Lennard-Jones  parameters  were  obtained  from  Refs  (24)  and  (25).  0^, 
the  reduced  collision  integral,  is  a  slowly  varying  function  of  the 
parameter  kT/c .  The  reduced  collision  integral  represents  the 
departure  from  rigid  sphere  behavior,  is  calculated  from  a  curve  fit 
of  data  tabulated  in  Appendix  B,  Table  B-2  of  Ref  (24). 
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The  viscosity  of  a  mixture  of  gases  is  calculated  using  the 


relation 


.  -V 

mix  /  M 


(3.20) 


where 


ij 


1 

yTS 


'  H  ' 

>  -1/2 

f  > 

1/2 

M  1 

1/*  • 

1  +  TT 

1  + 

1 

M 

M 

J  . 

J  J 

, 

Thermal  Conductivity.  The  Chapman- Enskog  I'^rmula  for  the  thermal 
conductivity  of  a  monatomic  gas  is 


k  - 

1 


1.9891  X  10 


Q 

1  k 


(3.21) 


where  0-0 
k  M 

The  calculation  of  the  thermal  conductivity  for  a  polyatomic  gas 
must  account  for  vibrational  and  rotational  excitation  as  well  as 
translational  energy.  The  Eucken  approximation  provides  a  method  to 
calculate  the  thermal  conductivity  of  a  polyatomic  gas. 

"  1 1  )  ^ 

The  formula  for  calculating  the  thermal  conductivity  of  the  mixture  is 
similar  to  Eq  (3.20) 
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k 

mix 


(3.23) 


-I 


X  k 
1  1 


t  xj> 


J-1 


ij 


Diffus ivltv.  The  expression  for  the  bimolecular  diffusivity  of 
species  i  into  species  j  is  very  similar  to  the  expressions  for  the 
viscosity  and  thermal  conductivity  of  of  a  species. 


(3.24) 


Note  that  the  Lennard- Jones  potential  is  now  the  intermolecular 
potential  between  one  molecule  of  species  i  and  one  molecule  of  species 

J. 


• 

12 

’  "I 

6 

-  4e 

iJ 

Q 

iJ 

.  , 

- 

a 

iJ 

the  reduced  collision  integral  for  collisions  between  unlike 
molecules,  is  calculated  from  a  curve  fit  of  data  tabulated  in  Appendix 
B,  Table  B-2,  (24).  The  Lennard-Jones  parameters  are  approximated  by 


a  -(a+o  )/2 

ij  i  J 


The  effective  diffusion  coefficient  of  species  i  into  the  mixture  of 
gases  is  obtained  by  assuming  that  the  diffusional  velocities  for  all  of 
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the  species  expect  1  are  equal  to  some  constant.  The  expression  for  the 
effective  diffusion  coefficient  is 


D  - 

i 


M 

I 

J-1 


(3.25) 


Chemical  Reactions 

Standard  AFWL  COIL  Rate  Package .  As  discussed  in  the  Introduction, 
AFWL  hosted  a  conference  to  review  the  current  state  of  the  Chemical 
Oxygen-Iodine  Laser  (COIL)  data  base  in  August  of  1987  (2).  One  product 
of  this  conference  was  a  new  standard  reaction  set  for  modeling  COIL 
devices .  The  standard  reaction  set  Includes  50  reactions ,  many  of  which 
appear  to  have  little  influence  at  standard  COIL  operating  conditions. 
As  part  of  the  preparation  for  the  kinetics  conference,  Captains  Glen 
Perram  and  Robert  Crannage  performed  a  set  of  sensitivity  analyses  to 
identify  a  reduced  set  of  equations  for  use  in  numerical  modeling. 
Their  work  identified  two  reduced  sets.  The  first  includes  22  reactions 
and  gives  excellent  agreement  with  the  full  reaction  set.  The  second 
includes  the  11  most  critical  reactions  and  gives  good  agreement.  In 
the  interest  of  computational  efficiency,  the  11  reaction  set  was  used 
in  this  investigation.  Table  1  lists  the  reactions  along  with  the 
accepted  empirical  rate  constants.  Note  that  this  is  the  same  set  of 
reactions  outlined  in  Chapter  II. 

There  are  several  points  that  should  be  made  about  the  reactions 
listed  in  Table  1.  First,  note  that  each  reaction  is  presented  in  terms 
of  a  forward  rate  constant  only.  Where  the  backward  reaction  is 
significant,  it  is  listed  as  a  separate  reaction  with  a  corresponding 
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TABLE  I 


Air  Force  Weapons  Laboratory  Reduced  COIL  Rate  Package 

11  Reaction  Set 


Kinetic  Process  Recommended  Rate 

Coefficient 
(cm  /molecule -sec) 


1. 

0,(^A) 

+ 

o/'a)  - 

- » 

02(^2) 

+ 

0 

V 

M 

k 

1 

- 

2.7x10*” 

2. 

0,(^2) 

+ 

HO  - 

- » 

0,(^2) 

+ 

H  0 

k 

6.7x10'” 

2 

2 

2 

2 

2 

3. 

I,(X) 

+ 

OjC'x)  - 

I  +  I 

+ 

0,(^S) 

k 

3 

- 

4.0x10*” 

4. 

I,(X) 

+ 

0  (^A)  - 

- » 

I* 

+ 

0,(^2) 

k 

7.0x10*” 

2 

2 

2 

4 

5. 

I,(X) 

+ 

* 

1  - 

- ) 

I* 

+ 

I 

k 

3.8x10*” 

2 

2 

5 

6. 

I  * 

2 

+ 

0/^A)  - 

- > 

I  +  I 

+ 

0,('2) 

k 

B 

- 

3.0x10*” 

7. 

I  * 

+ 

HO  - 

_ ^ 

I 

+ 

H  0 

k 

3.0x10*” 

2 

2 

2 

2 

7 

8. 

I 

+ 

0/^A)  - 

- » 

* 

I 

+ 

0/'2) 

k 

8 

- 

7.8x10*” 

9. 

* 

I 

+ 

0  (^2)  - 
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I 

+ 

0,(^A) 

k 

. 

2.7x10*” 

2 

g 

10. 

* 

I 

+ 

0  (^A)  - 

- > 

I 

+ 

0,(^2) 

k 

1.1x10*” 

2 

10 

11. 

* 

I 

+ 

HO  - 

2 

- » 

I 

+ 

H  0 

2 

k 

11 

- 

2.0x10*” 
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rate  constant.  For  example,  reaction  7  is  the  laser  pumping  reaction 
and  reaction  8  is  the  reverse.  Second,  the  rate  listed  for  reaction  2 
is  not  a  physically  measured  rate.  For  the  system  of  excited  oxygen  and 
water  vapor,  the  known  rates  are  for  reaction  1  and  for  the  total 
removal  of  0^(V)  to  all  products.  In  this  reduced  set,  the  rate  for 
reaction  2  has  been  set  to  produce  the  same  overall  removal  rate. 
Finally,  a  limited  data  base  does  exist  for  the  temperature  dependence 
of  these  rates.  This  data  is  presented  in  Table  2. 

Chemical  Rate  Equations .  Each  of  the  reactions  in  Table  1  can  be 
written  in  the  general  form  (26) 

i/'X  +  i/'X  + . . ¥  +  t/”X  +  .  (3.25) 

1122  1122 


or  as 


1 


I 


i/'X 

■  • 


(3.26) 


■•1  s-l 

where  X  is  any  chemical  species,  i/'  is  the  stoichiometric  coefficient 

for  the  reactant,  and  i/”  is  the  stoichiometric  coefficient  for  the 

1 

product.  Note  that  in  this  notation,  X  appears  on  both  sides  of  the 
reaction. 

The  molar  rate  of  formation  of  any  one  of  the  products  can  be 
expressed  as 


d 

dt 


X 

i/"k  N  r  X  1  Y  X  ] 

•  cA^-  1-*  2 


z 

2 


(3.27) 


where  k  [-]  (moles/cm^sec)  is  the  molar  rate  constant  and  the  product 

C 

on  the  right  hand  side  is  taken  over  all  of  the  reactants.  For  single 
step  reactions  or  sets  of  elementary  processes,  the  exponents,  z^'s,  are 
equal  to  the  stoichiometric  coefficients.  Similarly,  the  disappearance 


33 


TABLE  II 

Available  Temperature  Dependent  Rate  Data  for  the  COIL  Kinetics 


1  n  e  ia*28  (700/T) 

k  -  9.5x10  e 

1 

1  1  /  (1660/T) 

k  -  1.4x10  e 

5 

k  -  2.33x10*®/  T 

6 

*8 

,  3.1x10  (-*03^) 

k  -  - = -  e 

9  T 


,  4.0x10  ^  (700/1) 

k  -  - ;; -  e 

10  „3 . 8 
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of  each  of  the  reactants  is 


fr[X  ]--i/kNrx  iHx  1 

dtl-gj  scA‘-iJI-2-' 


(3.28) 


The  total  molar  rate  of  creation  of  each  species  is  the  summation  of  the 
creation  and  disappearance  terms  for  that  species  summed  over  all  of  the 
reactions,  N 


ii  [  *.  ]  -  i  ■  ■'(O'ac  ^  ]  ‘[  xj 


(3.29) 


n«t  n “  1 


The  mass  rate  of  creation  of  species  j  is  then  simply 


W  -  W  Ir  [  X  ] 
j  j  dt  s  J 


(3.30) 
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IV.  Niiniftrlcal  Analysis 


Integration  Scheme 

Numerical  methods  for  integrating  the  boundary- layer  equations  have 
been  studied  extensively  over  the  past  four  decades.  Anderson  et  al. 
(27)  provides  an  excellent  review  of  many  of  the  most  common  methods. 
Selection  of  a  particular  scheme,  an  explicit  or  an  implicit  method, 
requires  careful  consideration  of  such  factors  as  the  nature  of  the 
problem,  efficiency  and  accuracy  of  the  scheme,  and  the  complexity  of 
the  method.  For  this  investigation,  a  simple  explicit  marching  scheme 
was  selected.  It  has  the  advantage  of  being  very  straight  forward  and 
adapts  well  to  reactive  flows.  The  major  draw  back  of  an  explicit 
scheme  is  the  limitation  required  on  the  marching  step  size.  While  an 
implicit  scheme  would  be  unconditionally  stable  and  usually  more 
efficient,  it  would  be  more  difficult  to  implement.  In  addition,  the 
species  creation  rates  adds  a  further  constraint  on  the  marching  step 
size.  This  term  would  tend  to  increase  the  number  of  iterations 
required  to  converge  each  streamwise  step  in  an  implicit  method,  thus 
decreasing  the  efficiency  advantage. 

Finite -Difference  Formulas 

The  finite-difference  approximations  to  the  partial  derivatives  of 
an  arbitrary  function,  f(x,'6),  are  expressed  here  as  truncated  Taylor 
series  of  f(x,\fr)  on  a  nonuniform  grid.  Figure  7.  shows  the 

configuration  and  nomenclature  of  the  finite-difference  cell  in  the 
stream,  (a),  at  the  centerline,  (b) ,  and  at  the  freestream  boundary, 
(c). 

Interior  Nodes .  The  partial  of  f(x,^)  with  respect  to  x  is 
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Figure  7.  Finite-Difference  Grid  Notation,  a)  Interior  node, 
b)  Nozzle  centerline,  c)  Outer  streamline. 
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evaluated  using  a  first  order  forward-difference. 


(4,1) 

The  first  and  second  partials  with  respect  to  are  evaluated  using 
second  order  central-differences. 


a 


fCx.V*) 


(A0  A0^  +  t\j>^ 


+  (4.2) 


^2  f(x»  -  2 


A\&_f(x,^A^^>^)  -  (A^  +A\&^)f(x,V')  +  A^^f (x,^-A^_) 

-  +  -  + 


+  0(AV>^)  (4.3) 


Boundary  Nodes .  The  problem  of  interest  here  is  a  two  dimensional 
parallel  jet.  For  these  initial  numerical  investigations,  it  is  further 
assumed  that  the  jet  is  a  free  jet,  that  is  a  jet  in  an  unconfined  free 
stream.  The  jet  is  symmetric  about  the  centerline  which  implies  that 
there  are  no  transverse  gradients  at  the  centerline.  If  the  centerline 
is  selected  as  the  origin  for 

f(x,0)  -  0  (4.4) 

Approximating  3f/3\&  with  a  central-difference,  Eq  (4.2),  about  the 
centerline  (see  Figure  7  (b)),  Eq  (4.4)  is  satisfied  for  A^  -  A0^  and 
f(x,\fr-A\t_)  -  f(x,^A\fr^).  When  this  result  is  substituted  into  the 
central-difference  for  3  f/3x  ,  an  expression  for  the  second  partial 
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with  respect  to  i>  at  the  centerline  is  obtained. 


^2  f(x.O)  -  2 


f(x,AV>^)  -  f(x,0) 


(4.5) 


As  stated,  the  freestream  for  a  free  jet  is  unconfined.  If  a 
streamline  is  selected  (line  of  constant  ij>)  sufficiently  far  away  from 
the  jet,  one  would  expect  the  gradients  transverse  to  the  streamline  to 
be  zero. 


(4.6) 


A  similar  expression  to  Eq  (4.5)  can  be  derived  at  the  outer 
streamline,  ^  ^  (see  Figure  7  (c)). 

ft 


5^2  -  2 


f(x,0  -AV-  )  -  t(x,V',  ) 

_ ft  _ ft 


(4.7) 


Explicit  Finite -Difference  Equations 

The  finite-difference  formulas  for  di/di>  and  9^f/90^  defined  above 
can  be  rewritten  in  terms  of  general  difference  operators. 

A  f(x,V>  )  +  B  f(x.V-  )  +  C  f(x.V-  ) 

^  c  /  L\  I’  A  1  B  1  c 

^  -  - 5 -  (4.8) 


f(x,^) 


2 


A^f(x,^^)  +  B^f(x,\&^)  +  C^f(x,V>^) 
D 


(4.9) 


Details  of  the  values  for  the  coefficients  A,  B,  C,  A,  B,  C,  D, 

1  1  1  2  2  2 

0  ,  ^  ,  and  0  are  given  in  Appendix  B. 

Ac  C 
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The  thin- shear -layer  equations  in  Eqs  (3.13),  (3.14),  and  (3.15) 

were  approximated  using  the  finite-difference  formulas  defined  above. 
In  terms  of  the  finite-difference  operators,  the  governing  equations 
become , 


u(x+Ax)  -  u(x)  +  Ax 


1  dp 
pu  dx 


A  .  X  Au 


A^u 


A^ 


(4.10) 


W^(x+Ax)  -  W^(x)  +  Ax 


1  A  ^  A  .  _  X  AW  .  A^W 

-  W  +  -T-T  ipD  pu)  — rj  +  pa  pu  - J 

pu  j  A\6  "  A^  ^3  ^^2'' 


(4.11) 


T(x+Ax)  -  T(x)  +  Ax 


1  ^ 

+ 

Au  1 

pcp  dx 

Cp 

4^  J 

Ah  1 

1 

r  A 

(kpu) 

J 

Cp 

L  ^ 

Z  pW  h 
pucp  ^  j  j 

J 


A^T  ■' 


A^*  -> 


(4.12) 


Numerical  Stability 

One  disadvantage  of  explicit  finite-difference  formulations  is  that 
a  maximum  step  size  limitation  must  be  imposed  in  order  to  obtain  a 
stable  solution.  Determining  the  exact  criteria  for  the  step  size 
limitation  for  Eqs  (4.10),  (4.11),  and  (4.12)  is  difficult  due  to  the 
nonlinear  nature  of  the  equations.  Normally,  in  the  absence  of  chemical 
reactions,  stability  of  the  finite-difference  equations  will  be 
controlled  by  one  of  the  diffusive  terms  (viscous,  thermal,  or  mass). 
In  reacting  flows,  however,  an  additional  difficulty  is  encountered. 
The  formation  of  each  of  the  species  is  a  chemical  rate  process  with  a 
corresponding  time  scale.  If  the  time  scale  of  the  streamwise  marching. 
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(Ax/u) ,  is  larger  than  the  time  scale  of  the  chemical  rate  process,  the 
solution  of  the  chemical  kinetics  can  become  unstable;  however,  even 
"fast"  chemical  processes  tend  to  be  rapidly  modulated  by  the  ability  of 
the  reactant  species  to  diffuse  to  the  reaction  front. 

The  stability  of  a  finite-difference  equation  is  determined  by- 
examining  the  effect  of  introducing  small  errors  (from  numerical  round 
off  for  instance)  into  the  numerical  solution.  If  the  error  is 
amplified,  the  equations  are  unstable.  The  most  common  approach  used  to 
Investigate  stability  is  the  von  Neumann  error  analysis.  If  the 
numerical  solution  to  a  finite-difference  equation  is  expressed  as  the 
sum  of  the  exact  solution  plus  an  error,  it  can  be  shown  that  any 
numerical  error  must  also  satisfy  the  deference  equation.  The  error 
distribution  is  then  expressed  in  terms  of  a  Fourier  series  with  a  wave 
number  which  is  dependent  on  the  grid  spacing.  The  amplification  of  an 
arbitrary  Fourier  component  is  determined  by  substituting  into  the 
finite-difference  equation. 

In  order  to  examine  the  stability  of  the  finite-difference 
formulation  employed  here,  the  effect  of  the  diffusion  terms  on  the 
stability  of  the  mass  diffusion  equations  is  addressed  by  ignoring  the 
species  creation  term. 


dV  d_  \  dV' 

ax  “  av>  L  “  . 


(4.13) 


where  a  -  pDpxx.  Application  of  the  von  Neumann  stability  analysis  to  Eq 
(4.13)  can  be  simplified  by  recasting  the  finite-difference 
approximation  for  the  diffusion  term  in  a  slightly  different  form. 
Instead  of  expanding  the  diffusion  term  using  the  chain  rule,  the 
diffusion  term  is  evaluated  as  follows. 
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l•&)l . -1-s-IL 

I  +  Aif>J  j 


(4.14) 


where. 


I  •".■II 


t  +  a  1  f  W  -  W 
J-n  J  j-*-!  J 

2  J  [ 


(4.15) 


(  aw  1 1  f  “j-i 


w  -  w 

j  j-1 


(4.16) 


The  finite-difference  approximation  to  Eq  (4.13)  is  then 


a  (W  -  W  )  a  (W  -  W  ) 

W  -  W  +  2  Ax  i-  +  J 

i-t-i  1  +A\J  )  Aij>  (Ai>+Ail>  ) 

+  +  -  -  +  • 


(4.17) 


The  von  Neumann  stability  analysis  can  then  be  applied  to  Eq  (4.17)  by 
approximating  the  error  as 


,  .  ax  lk\fr 

€(x,i>)  -  e  e 


and  requiring  the  absolute  magnitude  of  the  amplification  be  less  than 
or  equal  to  one.  The  resulting  stability  criteria  for  Eq  (4.17)  is 


Alt’  Arp Alp  +AV>^) 

Ax  <  ^  .  - ^--7 - r  X  C 

2  ( A\b  a  +  Alp  a  ) 
j  +  l/2  +  i-l/2 


(4.18) 


The  stability  criteria  for  momentum  or  thermal  diffusion  is  identical 
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with  a  -  fipu  and  a  -  k/)u,  respectively.  In  the  absence  of  large 
velocity  or  temperature  gradients,  mass  diffusion  usually  controls  the 
stability. 

This  analysis  is  only  approximate  and  does  not  provide  an  absolute 
stability  criteria.  Coupling  between  the  equations,  the  presence  of 
artificially  large  gradients  in  the  initial  conditions,  and  the 
stability  of  the  chemical  rate  processes  are  not  addressed.  The 
constant  C  (0.0  <  C  <1.0)  in  Eq  (4.18)  allows  the  diffusion  limit  to  be 
further  reduced  in  order  to  obtain  a  stable  solution.  Experience  has 
shown  that  C  0.05  is  required  for  the  problem  under  investigation 
here.  In  addition,  for  the  initial  starting  conditions  where  the 
species  gradients  are  set  to  step  functions  at  the  jet  boundary,  C  is 
further  reduced  to  0.01  and  allowed  to  Increase  linearly  with  streamwise 
position  to  0.05. 


43 


V.  Results  and  Discussion 


Introduction 

The  COIL  kinetic  rate  package  was  combined  with  the 
thin- shear -layer  equations  and  a  simple  explicit  integration  scheme 
(Chapters  III  and  IV)  to  create  a  numerical  model  to  be  used  in  concert 
with  the  planned  two  dimensional  jet  experiment  to  investigate  I^ 
dissociation.  The  resulting  code,  COILFT,  is  listed  in  Appendix  C. 
This  section  presents  the  results  of  several  numerical  investigations 
conducted  to  validate  this  code  and  demonstrate  its  capabilities  as  an 
investigative  tool.  The  first  case  presented  is  a  comparison  with  the 
results  of  a  nonequilibrium  chemistry  code  to  validate  the  chemistry 
model  installed  in  COILFT.  The  next  series  of  cases  compares  a  set  of 
premixed  cases  (i.e.,  the  0^  and  I^  are  uniformly  mixed)  and  simple 
jet-mixed  cases  (i.e.,  a  secondary  jet  of  I^  is  injected  in  a  parallel 
sheet  into  the  0^  primary  flow)  in  order  to  investigate  the  influence  of 
"mixing"  on  the  dissociation  process.  The  last  set  of  cases  are 
intended  to  demonstrate  the  role  that  water  vapor  plays  in  the 
dissociation  process  for  both  the  premixed  and  jet-mixed  cases.  The 
intent  is  to  perform  some  simple,  well-defined  numerical  experiments 
(simulations)  that  will  aid  in  building  some  rational  feel  for  and  raise 
questions  about  common  perceptions  concerning  the  performance  of  COIL 
lasers  (c.f.,  preceding  chapters).  Finally,  although  this  work  is 
intended  to  stand  on  its  own,  the  jet  dimensions  and  velocity  range  were 
chosen  to  bear  some  resemblance  to  the  concurrent  experimental  program. 
In  this  regard,  then,  the  results  of  this  study  will  be  useful  in 
guiding  the  eventual  experimental  design  of  the  concurrent  effort  (c.f., 
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previous  chapters) . 


Validation  of  COIL  Chemistry  Package 

Figure  8  shows  a  comparison  of  the  time  dependent  concentration 
profiles  of  •  ^2’  using  COILFT  and  those  from 

the  presently  accepted  "canonical"  nonequilibrium  chemistry  code, 
SENCHEM  (28).  SENCHEM  calculates  the  concentrations  of  various  mixtures 
due  to  the  eleven  reaction  set  referred  to  earlier  (c.f.,  Table  I, 
Chapter  III)  as  a  function  of  time  at  constant  pressure  and  temperature. 
It  should  be  noted  the  SENCHEM  is  purely  a  kinetics  code  and  is 
incapable  of  describing  fluid  mechanic  coupling.  In  order  to  compare  to 
the  SENCHEM  results,  COILFT  was  set  up  with  a  uniformly  mixed  (i.e.  no 
mixing)  initial  condition  for  the  same  mixture,  pressure,  and 
temperature  as  SENCHEM  and  a  uniform  initial  velocity.  The  resulting 
concentration  of  each  species  verses  position  downstream  was  then 
transformed  to  time  by  dividing  the  position  by  the  velocity.  The 
initial  conditions  used  for  this  comparison  are  listed  in  Table  III. 
This  comparison  between  the  two  computational  models  was  used  to  trouble 
shoot  the  installation  of  the  COIL  chemistry  package  into  the 
thin- shear -layer  code  and  to  validate  the  results.  Reference  to 
Figure  8  shows  that  essentially  no  differences  exist  between  the  results 
of  the  two  codes,  i.e.,  we  can  have  high  confidence  that  the  kinetics  in 
COILFT  function  properly. 

Mixing  Studies 

The  next  six  cases  presented  here  are  intended  to  provide  insight 

into  the  role  mixing  plays  in  the  I^  dissociation  process.  The  first 

three  runs  were  for  parallel  jet  cases  with  varying  ratios  of  the 

secondary  jet  (I  and  diluent)  velocity  to  the  primary  freestream  (0 
2  2 
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SENCHEM  and  COILFT. 


table  III 


Initial  Conditions  for  COIL  Chemistry  Validation  Case 


Temperature 

Pressure 

Velocity 

300.0  K 

3 . 33x10*  dyne/cm^ 

5,550.0  cm/sec 

(25.0  torr) 

Concentrations 

Species 

Number  Density 
(molecule/cm  ) 

Mole  Fraction 

8.85x10^® 

0.1100 

0^<V) 

5.91xl0“ 

0.0734 

8.05xl0“ 

0.0001 

I,(x) 

3.22xl0“ 

0.0040 

I* 

0.0 

0.0 

2 

I 

0.0 

0.0 

* 

I 

0.0 

0.0 

He 

6.39x10^^ 

0.7936 

H  0 

3.22x10^’ 

0 . 0440 

Cl 

1.12xl0“ 

0.0149 

2 
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and  diluent)  velocity.  At  the  lower  limit,  this  ratio  is  set  to  one  and 
mixing  between  the  two  streams  is  due  solely  to  molecular  diffusion  of 
the  species.  As  the  jet  velocity  increases  relative  to  the  freestream 
velocity,  the  oxygen  stream  becomes  "entrained"  into  the  iodine  stream 
resulting  in  increased  mixing.  Entrainment  is  used  here  to  describe  the 
process  by  which  the  higher  velocity  secondary  jet  accelerates  the 
primary  flow  through  shear;  the  increasing  velocity  of  the  primary  flow 
gives  rise  to  a  component  of  velocity  in  the  primary  flow  toward  the 
secondary  jet,  thereby  adding  a  convective  transport  mechanism  for  the 
O^(^A)  in  the  direction  of  the  I^.  This  transport  tends  to  increase  the 
concentration  gradient  for  molecular  transport,  thus  enhancing  the  rate 
of  molecular  mixing.  The  next  three  are  "premixed"  cases  representing 
the  ideal  instantaneous  mixing  of  the  oxygen  and  iodine  streams  of  the 
three  jet  cases  at  the  initial  plane,  not  unlike  the  chemistry 
verification  run  of  the  previous  section. 

Geometry  and  Initial  Conditions.  The  geometry  and  flow  conditions 


for  the  jet-mixed  cases 

were 

chosen  to 

approximate  those 

of  the 

concurrent  experimental 

study 

discussed 

in  Chapter  II. 

The 

configuration  and  expected  performance  of  both  the  oxygen  generator  and 
iodine  injection  nozzle  are  discussed  briefly  Appendix  A.  Figure  9  is  a 
schematic  of  the  computational  domain  for  the  parallel  jet  cases 
(simplified  from  the  experimental  test  section  diagram  in  Figure  Al, 
Appendix  A)  which  illustrates  the  setup  of  the  numerical  cases.  The 
solutions  were  begun  at  the  nozzle  exit  plane  where  the  velocity, 
temperature,  and  species  mole  fractions  were  taken  to  be  uniform  "slug" 
profiles  for  both  the  jet  and  freestream  flows.  The  pressure  was 
specified.  Taking  advantage  of  the  symmetry,  the  numerical  solutions 
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Figure  9 , 


Schematic  of  the  Computational  Domain  for  the  Parallel  Jet 
Cases . 


were  calculated  from  the  nozzle  centerline  to  the  "wall"  (the  wall 
itself  was  not  modeled) .  The  combined  height  (cavity  half  height)  of 
the  jet  and  freestream  flows  was  0.625  inches  (1.59  cm)  and  the  nozzle 
height  for  the  secondary  jet  alone  was  0.0375  inches  (0.095  cm).  Six 
finite -difference  nodes  (in  the  y/\&  direction)  were  placed  in  jet  with 
twenty-seven  in  the  freestream. 

Three  jet-mixed  cases  with  ratios  between  the  jet  velocity  and  the 
freestream  velocity  of  1/1,  3/1,  and  6/1  were  run.  The  jet  Mach  Numbers 
ranged  from  0.1  to  0.4.  The  freestream  Mach  Number  was  approximately 
0.1.  The  velocity  ratio  of  the  jet  was  changed  by  changing  the  amount 
of  diluent  (helium  carrier  in  this  case)  in  the  jet.  Note,  the  flow 
rate  of  iodine  was  held  constant  although  the  the  iodine  mole  fraction 
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changed  with  increasing  diluent.  The  freestream  flows  were  held 
constant.  The  initial  conditions  for  these  three  cases  are  listed  in 
Table  IV. 

For  comparison,  three  premixed  cases  were  also  run.  These  cases 
were  created  by  uniformly  mixing  the  freestream  and  jet  flows  for  each 
of  the  jet-mixed  cases.  In  order  to  allow  for  some  rational  basis  of 
comparison  between  the  premixed  and  jet-mixed  cases,  the  velocities  for 
the  premixed  cases  were  set  equal  to  the  corresponding  jet  velocity.  By 
comparing  premixed  and  jet -mixed  cases  in  this  way,  the  common  factor 
was  the  time  that  the  was  exposed  to  the  O^(^A).  The  criteria  was 
not  strictly  adhered  to  however,  since  the  secondary  jet  did  lose 
velocity  with  time.  The  initial  conditions  for  these  premixed  cases  are 
listed  in  Table  V. 

Assumptions  and  Simplifications .  Before  examining  the  results  of 
these  test  cases,  it  is  important  to  be  aware  of  the  assumptions  and 
simplifications  made  in  selecting  both  the  initial  conditions  and 
boundary  conditions  for  these  numerical  cases.  Some  of  the  caveats 
placed  on  these  results  are  implicit  in  the  solution  while  others  have 
been  made  in  order  to  isolate  the  specific  influence  of  mixing  and  water 
vapor  and  to  make  the  results  easier  to  interpret.  Where  prudent,  the 
boundary  conditions  have  been  simplified  to  reduce  the  computational 
complexity. 

An  implicit  assumption  in  the  discussion  to  this  point  is  that  the 
chemical  kinetics  package  employed  in  COILFT  accurately  reflects  the 
actual  physical  processes.  Confirming  the  accuracy  of  the  current  COIL 
kinetics  package  is,  in  fact,  the  purpose  of  the  synergistic  concurrent 
experimental  and  computational  program  of  which  this  effort  is  a  part. 


50 


Initial  Conditions  for  the  Jet  Mixing  Cases 
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TABLE  V 


Initial  Conditions  for  the  1-D,  Premixed  Cases 


Associated 

Jet-Mixed  Case 

1/1 

3/1 

6/1 

Temperature  (K) 

400.0 

400.0 

400.0 

2 

Pressure  (dyne/cm  ) 

2.9795x10* 

3.1449x10* 

3.3635x10* 

Velocity  (cm/sec) 

6228.4 

j 

17702.5 

31273.8 

Concentrations 
(mole  fraction) 

Species 

O^('Z) 

0.1435 

0.1322 

0.1195 

Oj('a) 

0.0957 

0.0881 

0.0797 

OjC'S) 

0.0 

0.0 

0.0 

Iju) 

0.0048 

0 . 0044 

0 . 0040 

I  ♦ 

2 

0.0 

! 

0.0 

0.0 

I 

0.0 

0.0 

0.0 

* 

I 

0.0 

0.0 

0.0 

He 

0.7560 

0.7753 

0.7968 

H  0 

2 

0.0 

0.0 

0.0 

Cl 

2 

0.0 

0.0 

0.0 
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Even  if  the  individual  kinetic  rates  are  correct,  the  reduced  eleven 

reaction  set  is  not  necessarily  valid.  Note,  for  instance,  that  the 

initial  concentrations  of  H  0  and  0  (^2)  are  set  to  zero  for  all  of  the 

2  2 

cases  in  Table  IV  and  V.  In  an  actual  device,  it  is  essentially 
impossible  to  eliminate  all  of  the  H^O  vapor.  Even  at  low 

concentrations,  H^O  vapor  is  the  strongest  quenching  mechanism  for 

0  (^E),  I  and  I  .  This  fact  allowed  the  removal  of  all  other 

A  2 

quenching  mechanisms  from  the  eleven  reaction  set.  In  the  absence  of 

H  0,  other  loss  mechanisms  not  included  in  Table  I  could  become 
2 

important.  Furthermore,  assuming  that  our  knowledge  of  the  complete 

reaction  set  is  sound,  we  know  that  in  the  presence  of  H^O,  O^(^E) 

reaches  an  equilibrium  value  between  production  from  O^(^A)  pooling  (Eq 

2.3)  and  quenching  with  HO  (Eq  2.4).  In  the  absence  of  HO,  the  0  (^E) 

2  2  2 

concentration  would  slowly  rise  to  unrealistic  concentrations.  Since 
O^(^E)  plays  a  role  in  initiating  the  dissociation  of  I^  (Eq  2.6),  the 
dissociation  rates  would  be  accelerated.  Ignoring  H^O  in  the  initial 
conditions  effectively  slows  this  affect.  The  assumption  of  no  H^O  in 
the  initial  conditions  does,  however,  allow  for  a  clearer  interpretation 
and  makes  understanding  of  the  results  presented  here  possible.  The 
role  of  H^O  is  treated  independently  in  a  separate  set  of  numerical 
cases . 

The  most  important  assumption  made  here  is  that  the  jet  flow  is 
laminar.  The  available  transition  data  base  is  sparse.  Finding  clear 
transition  data  for  this  problem  is  further  complicated  by  the  presence 
of  temperature  and  concentration  gradients,  chemical  reactions,  and,  in 
the  case  of  the  concurrent  experiment,  a  confining  wall.  The  Reynolds 
number  based  on  the  jet  height,  R  ,  for  the  6/1  jet  case  is 

d 
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approximately  200.  Birkhoff  et.  al.  (16)  suggest  that  for  R  <  1000 

d 

(round,  free  jets),  the  jet  should  remain  laminar  (molecular  transport 
only),  however,  for  R  >  100  the  development  of  periodic  large  scale 
structure  is  certainly  a  possibility.  The  whole  issue  is  further 
complicated  by  the  use  of  the  words  laminar  and  turbulent,  where  in  many 
cases  large  scale  vortical  structure  are  not  considered  turbulent 
because  they  have  not  cascaded  down  to  the  higher  frequency,  small  scale 
eddies.  The  assumption  of  laminar  flow,  however,  is  pragmatic  from  a 
computational  point  of  view  and  may,  in  fact,  be  correct.  Ultimately, 
the  validity  of  this  key  assumption  must  be  confirmed  or  rejected  based 
on  experimental  evidence . 

It  should  also  be  stressed  that  this  study  represents  a  set  of 
"numerical  experiments,"  and  as  such,  is  not  meant  to  simulate  an  actual 
experiment  in  anymore  than  a  general  way.  This  is  apparent  from  the 
initial  conditions.  Reference  to  Figure  Al,  Appendix  A,  clearly 
demonstrates  that  an  actual  experiment  would  involve  the  growth  of 
boundary  layers  on  the  inner  and  outer  walls  of  the  nozzle.  The  no-slip 
condition  requires  that  the  velocity  at  the  edge  of  the  nozzle  blade 
(the  Interface  between  the  two  streams)  be  zero  although  concentration 
gradients  would  be  large.  In  the  numerical  model  the  initially  large 
gradients  between  the  freestream  and  jet  "slug"  profiles  are  smeared  out 
over  the  two  adjacent  grid  nodes  between  the  jet  and  the  freestream. 

One  would  also  expect  the  growth  of  a  boundary  layer  on  the  cavity 
walls.  This  viscous  shear  loss  coupled  with  heat  release  due  to 
chemical  reaction  should  cause  some  pressure  gradient  in  the  st..eamwise 
direction.  Since  the  freestream  Mach  number  is  small  (M  »  0.1)  and  the 
amount  of  heat  released  is  also  small  (near  resonant  energy  transfers), 
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the  pressure  gradient  would  be  due  to  the  shear  on  the  wall.  As 
mentioned  in  Chapter  4,  no  attempt  has  been  made  in  this  numerical  study 
to  model  the  streamwise  pressure  gradient.  The  pressure  gradient 
instead  was  set  to  zero  for  all  of  the  cases.  Because  the  solution  is 
being  calculated  along  lines  across  which  no  net  mass  flows,  i.e., 
constant  the  error  due  to  this  assumption  can  be  qualitatively 
evaluated  by  examining  the  variation  in  the  y- coordinate  of  the  outer 
streamline.  The  maximum  change  was  8%  for  the  6/1  jet-mixed  case.  If 
required,  a  streamwise  pressure  gradient  could  be  specified  based  on 
measured  wall  pressure  data  (See  Figure  Al,  Appendix  A).  The 
boundary  layer  at  the  wall  could  also  be  modeled  by  imposing  a  no-slip 
condition. 

General  Results .  The  output  of  the  numerical  code,  COILFT, 
provides  information  on  both  the  flow  field  properties  and  species 
concentrations.  Figures  10  and  11  show  the  development  of  the 
transverse  velocity  and  temperature  profiles  for  the  6/1  velocity  ratio 
Jet  case  as  a  function  of  xp.  The  slug  profile  initial  conditions  at  the 
nozzle  exit  plane  quickly  develop  into  a  series  of  profiles  which  may  be 
compared  with  a  self-similar  free  jets  (19).  Although  the  comparison  is 
complicated  by  dissimilar  species,  chemical  diffusion,  and  reactions, 
the  evolution  of  the  velocity  and  temperature  profiles  closely  mimic 
those  of  the  self-similar  jet  (19),  and  give  high  confidence  that  the 
code  performed  properly  with  regard  to  simulating  the  fluid-mechanic 
environment.  Figure  12  shows  the  development  of  the  species  profiles 
transverse  to  the  flow,  three  centimeters  downstream  of  the  nozzle  exit 
plane.  Note  that  the  freestream  profiles  are  undisturbed  confirming 
that  ignoring  the  Influence  of  the  cavity  walls  was  a  reasonable 
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Figure  11.  Development  of  the  Transverse  Temperature  Profiles  as  a  Function 
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Figure  12.  Transverse  Species  Concentration  Profiles  for  the  6/1  Jet  Case 
3  cm  Downstream  from  the  Nozzle  Exit  Plane. 


assumption.  Figure  13  shows  the  development  along  the  nozzle  centerline 

in  the  streamwise  direction.  Close  inspection  of  Figures  12  and  13 

★ 

shows  the  diffusion  of  0^  into  the  jet  and  the  rise  of  both  I  and  I  as 
the  dissociation  process  proceeds.  Although  these  detailed  maps  of  the 
individual  species  concentrations  will  become  important  when  attempting 
to  compare  numerical  results  with  experimental  data,  developing  a  "feel" 
for  the  dissociation  process  from  these  profiles  alone  is  often 
difficult  because  it  requires  examining  the  changes  in  several  different 
species  as  a  function  of  both  the  streamwise  and  transverse  position  at 
the  same  time.  To  aid  in  developing  a  "feel"  for  the  process  the  flow 
field  properties  and  the  species  concentrations  were  manipulated  to  form 
Integrated  parameters  which  are  more  characteristic  of  the  overall 
product  of  the  dissociation  process. 

StTPArawi RP  Comparisons .  The  two  questions  of  most  concern  to  this 
investigation  were  what  is  the  rate  of  dissociation  and  how 

efficiently  does  the  dissociation  occur.  Integration  of  the  transverse 
profiles  of  both  the  0^  and  as  the  flow  develops  in  the  streamwise 
direction  provided  a  means  of  addressing  both  of  these  questions. 

The  fraction  of  the  I  dissociated  can  be  defined  in  terms  of  the 

2 

total  mass  flow  rates  of  iodine  through  a  plane  perpendicular  to  the 
nozzle  centerline  as 


Am(I  )  ni(I  ) 

Diss - = - 1 - — 


m(I  ) 

2  TOT 


md  ) 

2  TOT 


(5.1) 


where 


m(I  )  -  m(I  )  +  m(I 
2  2 


2 
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Figure  13.  Streamwise  Species  Concentration  Profiles  for  the  6/1 
the  Nozzle  Centerline. 


m(IJ _ -  m(I„)  +  +  m(I)  +  n>(I*) 


2  TOT 


Total  dissociation  corresponds  to  Diss  -  1.  The  total  mass  flow  rate  of 
species  i  through  a  plane  was  found  by  integrating  the  product  of  the 
mass  flow  rate  (,r/>)  and  the  mass  fraction  of  species  i  across  the  flow. 
This  integral  was  approximated  on  the  finite -difference  grid  as 


"  (  ) 


A0 


(5.2) 


J-i 


where  n  is  the  total  number  of  nodes  transverse  to  the  flow. 

Figure  14  shows  the  dissociation  fraction  for  the  three 

jet-mixed  cases  and  the  corresponding  premixed  cases  as  a  function  of 

the  streamwise  position.  Note  that  initially,  just  as  for  the  RDA 

results,  Figure  4,  Chapter  II,  the  rate  of  dissociation  for  the  premixed 

cases  is  slower  than  for  the  jet-mixed  cases.  Comparison  of  the  species 

concentrations  as  a  function  of  the  distance  downstream  for  the  premixed 

case  in  Figure  15  with  those  at  the  jet  boundary  for  the  6/1  jet-mixing 

case  in  Figure  16  shows  a  higher  initial  I  concentration  with  a 

* 

corresponding  rapid  rise  in  the  concentrations  of  I ,  I  ,  and,  most 
importantly,  I^^ .  The  result  is  a  much  more  rapid  rate  of  dissociation 
at  the  jet  boundary.  The  drop  off  in  I^^  concentration  near  the  2  cm 
point  occurs  simply  because  the  I^  is  being  depleted.  Because  in  the 
premixed  case,  all  species  are  completely  mixed,  the  initial 
concentration  of  I^  is  lower  for  the  premixed  case.  The  rate  of 
formation  of  all  of  the  excited  species  in  the  premixed  case  shows  a 
much  more  gradual  rise.  As  the  I  concentration  builds,  however,  the 
rate  of  dissociation  increases  due  to  a  hand  off  between  the  slower 
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dissociation  mechanisms,  Eqs  (2.6)  or  (2.12),  and  (2.14),  and  the  faster 
mechanism  of  Eqs  (2.13)  and  (2.14)  as  discussed  in  Chapter  II.  All  of 
the  jet-mixed  cases  display  a  slow  roll  off  in  the  dissociation  rate  and 
the  premixed  cases  eventually  overtake  the  jet  cases.  Note,  that  the 
6/1  case  which  corresponds  to  highest  entrainment  rate  stays  ahead  of 
the  premixed  case  for  the  longest  period  of  time. 

In  Chapter  I,  the  efficiency  of  the  I^  dissociation  process  was 
defined  in  terms  of  the  parameter  N,  the  number  of  O^(^A)  molecules 
required  to  dissociate  one  molecule.  An  efficiency  parameter  similar 
to  N  can  be  defined  for  the  numerical  results  in  terms  of  the  total 
(integrated  across  the  flow)  loss  of  the  excited  "chemical  energy"  in 
the  flow  divided  by  the  change  in  the  flow  rate  of  I^.  Initially,  all 
of  the  chemical  energy  in  the  flow  is  in  0^(V). 

E  -  Ef  0  Oa)  ) 
in  '  2  •'in 

E  is  the  total  energy  flux  [-]  (kJ/sec)  in  the  flow,  n  is  the  total 

molar  flow  rate  [-]  (mole/sec),  and  e  is  the  energy  per  mole  [-] 

(kJ/mole).  e  corresponds  to  the  excitation  energy  in  Eq  (2.18),  Chapter 

II  (also,  see  Figure  1).  The  subscript  in  refers  to  the  initial  plane. 

At  any  streamwise  location,  the  energy  still  available  in  the  flow  is 

1  * 

defined  as  the  energy  in  both  0^(  A)  and  I  .  Any  energy  that  transfers 

to  0  (^2)  or  I  ^  is  considered  as  a  loss. 

2  2 
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E  -  £(  O^(^A)  )  +  eC  I*  ] 

-  A[  0^(V)  )  «[  O^(^A)  )  +  n[  I*  )  «[  I*  3  (5.4) 

The  total  energy  loss  between  the  initial  plane  and  any  streamwise 
location  is  then 


AE  -  E  -  E 

IN 


which  can  be  expressed  in  terms  of  the  transition  energy  of  O^(^A) 


AE 


^[1*3 


«(  0,(  A)  3 


—  -  n[  o/a)  3^^  -  nC  O^(^A)  3  -  - - ^  n(  I*  3  (5.5) 


0,(  A)  3 


The  total  loss  of  to  dissociation  can  be  computed  by  accounting  for 

both  I  and  I  ^  as 
2  2 

An(  3  -  -  nC  ^2  )  -  n[  ^ 


The  efficiency  can  then  be  defined  as 

ae/«C  0  ('a)  3 

N(x)  -  — L -  (5  7) 

Ar,(  ) 

Note,  that  the  lower  this  number  is,  the  more  efficient  the  dissociation 
process  is  (less  energy  lost).  The  theoretical  limit  based  on  energy 
required  to  dissociate  is  2  (See  the  energy  level  diagram  in 

Figure  1,  Chapter  I).  For  initially  very  small  dissociation  fractions, 
the  denominator  of  Eq  (5.7)  is  near  zero,  which  causes  erratic  behavior 
in  computed  efficiencies  from  the  code;  these,  however,  settle  out 
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rapidly  as  An(l^  becomes  sufficiently  large  to  overcome  the  roundoff 

errors  involved  in  subtracting  nearly  equal  nvunbers. 

Figure  17  shows  the  streamwise  variation  of  the  efficiency  factor 

defined  in  Eq  (5,7)  for  the  six  cases  with  the  curves  beginning  after 

the  calculated  error  has  settled  out.  The  slight  difference  in 

efficiency  between  the  premixed  and  jet-mixed  cases  can  be  seen  in 

Figure  17.  The  slight  increase  in  N(x)  (i.e.,  decrease  in  efficiency, 

c.f.  above)  in  the  streamwise  direction  is  due  to  the  relatively  slow 
1  * 

0^(  A)  and  I  pooling  reactions,  Eqs  (2.2)  and  (2.4).  Keep  in  mind  that 

with  no  HO  in  the  initial  conditions,  there  are  no  other  loss 
2 

mechanisms,  so  that  the  value  of  the  efficiency  factor  should  be  close 
to  the  minimum  limit  of  2.  Although  defined  similarly,  the  efficiency 
factor  from  Eq  (5.7)  illustrates  an  important  point  about  the  N  commonly 
used  in  Eq  (1.1);  N  is  defined  as  the  number  of  0^(V)'s  required  to 
dissociate  one  at  100%  dissociation.  N  is  then  an  asymptotic  limit 
of  the  factor,  N(x)  in  Figure  17  as  the  dissociation  proceeds  to 
completion.  Difficulty  arises,  however,  in  defining  N  precisely. 
Examination  of  Eq  (5.7)  and  Figure  17  shows  that  even  after  the 
dissociation  has  reached  100%  the  loss  mechanisms  will  continue  to  cause 
the  apparent  N  to  increase.  This  example  points  out  that  N,  as  defined 
and  used  in  Eq  (1.1),  is  not  a  universal  constant  for  COIL  devices  but 
tied  to  a  specific  device  and,  in  fact,  to  a  specific  run  condition. 
Clearly,  N  will  be  of  little  value  when  making  large  excursions  from  the 
current  baseline  designs  and  operating  conditions. 

Transverse  Comparisons .  For  thin  shear  layers,  gradients 
transverse  to  the  flow  direction  are  more  important  to  the  development 
of  the  flow  than  gradients  in  the  streamwise  direction,  since  the 
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transverse  profiles  provide  an  indication  of  the  structure  of  the  jet  . 
Once  again,  combined  parameters  which  characterize  COIL  system 
performance  provide  a  clearer  picture  than  attempting  to  look  at  the 
individual  profiles.  The  first  parameter  of  interest  is  the  local 
dissociation  fraction.  The  local  dissociation  fraction  is  the  mass 
fraction  of  atomic  iodine  at  any  point  in  the  flow  divided  by  the  total 
mass  fraction  of  iodine  (atomic  or  molecular)  at  the  point. 


I 

2 


Diss  - 


W(  I  ) 

W[  I  ) 

Z  ■’  TOT 


where  W  is  the  mass  fraction  and 


(5.8) 


wC  I  -  WC  I  ]  +  W(  I*  ) 


(5.9) 


and 


^2  \oT  -  W(  Ij  +  W[  I*  )  +  W(  I  }  +  w[  I*  }  (5.10) 

Complete  dissociation  corresponds  to  Diss  -  1.  As  before,  some 

confusion  arises  because  the  I  dissociation  is  not  defined  where  the 

2 

total  mass  fraction  of  iodine  is  zero  or  well  behaved  where  it  is  small. 
Plotting  the  total  mass  fraction  of  aids  in  interpreting  the 
dissociation  fraction. 

The  O^(^A)  "yield",  fraction  of  0^  that  is  O^(^A),  along  with  the 

total  mass  fraction  of  0^  provides  means  of  tracing  the  "entrainment"  of 

the  0^  from  the  primary  freestream.  The  yield  (c.f.,  above)  is  defined 

as  the  fraction  of  0  (^A)  in  the  total  0 

2  2 
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YIELD  - 


(5.11) 


W(  O^('a)  ) 

W(  0  ] 

V  2  ^  TOT 

where 

W(  0^  -  W(  O^('s)  )  +  w(  O^('a)  )  +  W(  O^('e)  )  (5.12) 

Finally,  the  "gain"  provides  an  indication  of  where  the  atomic 
iodine  has  reached  an  inverted  population.  As  discussed  in  Chapter  II, 
Eq  (2.7),  inversion  occurs  when 

[I*] 

-  >  0.5  (5.13) 

[I] 

The  gain  is  defined  as 

g  -  o  (  [  I*  ]  -  0.5  [  I  ]  )  (5.14) 

where  a  is  the  stimulated  emission  cross  section,  a  -  5.7X10  cm  . 

Figure  18  shows  the  transverse  profiles  for  the  total  I^  mass 

fraction,  Eq  (5.10),  local  I^  dissociation  fraction,  Eq  (5.8),  total  0^ 

mass  fraction,  Eq  (5.12),  0^(V)  traction,  Eq  (5.11),  and  gain, 

Eq  (5.14),  at  a  streamwise  position  1  cm  from  the  nozzle  exit  plane  for 

the  6/1  jet  case.  The  origin  for  the  y  cooriinate  is  the  nozzle 

centerline,  and  the  initial  boundary  between  the  et  and  freestream 

flows  is  at  y  -  0.09  cm;  the  profiles  are  displayed  out  to  y  -  0.8  cm, 

approximately  half  of  the  computational  domain.  The  remaining 

freestream  flow  is  undisturbed  by  the  Jet  expansion.  Figure  18  shows 

that  the  0^  and  I^  streams  have  only  just  began  to  mix  by  the  1  cm 

location.  The  I  dissociation  fraction  has  been  deliberatelv  set  to 
2 


70 


2.0 


3  ^  Q 


> 


c 


’(C 


CD 


(uuo/o^)  ‘ureo 


cf 


o 


(0 


o 

a 


P|9!A  fO)M  ‘ssiQ  ®l  f DM 


71 


Figure  18.  Transverse  0^  Mixing,  Dissociation,  Gain  Profiles  for  the  6/1 
Jet  Case  at  x  -  1  cm. 


zero  for  total  mass  fractions  less  than  0.01  to  avoid  any  confusion 
from  artificially  high  dissociation  fractions  in  regions  where  the 
iodine  concentrations  are  numerically  very  small.  This  accounts  for  the 
dramatic  drop  off  in  the  dissociation  fraction  near  y  -  2.5  cm.  The 
is  only  slightly  dissociated  in  the  mixing  region  (  »  0.15  ).  No 
dissociation  has  occurred  in  the  core  of  the  jet  by  the  this  point. 
Note,  gain  occurs,  according  to  the  arguments  in  Chapter  2  and  by 
examination  of  the  definition  in  Eq  (2.14),  only  where  has 

dissociated  to  create  I  atoms  and  where  the  O^(^A)  fraction,  Eq  (5.11) 
is  greater  than  0.15.  The  gain  profile  in  Figure  18  illustrates  this 
point . 

Figures  19  and  20  repeat  the  profiles  of  Figure  18  for  the  6/1  jet 
case  at  the  x  -  3  cm  and  x  -  6  cm  planes  respectively.  As  the  flow 
progresses  downstream,  the  in  the  jet  diffuses  outward  while  the  0^ 
in  the  jet  diffuses  inward  aided  by  convective  transport  as  the  0^  is 
"entrained"  into  the  faster  moving  jet  flow.  Although  the  0^ 
concentration  on  the  jet  centerline  increases  rapidly,  the  O^(^A) 
fraction  drops  quickly  to  zero  in  the  jet  core.  The  which  diffuses 
into  the  O^(^A)  rich  freestream  is  dissociated  rapidly.  The  largest 
concentration  of  I^,  however,  remains  in  the  jet  core  even  at  x  -  6  cm 
and  nearly  one  half  of  the  I  in  the  core  is  still  undissociated.  Note 
the  location  of  the  gain  regions.  Once  again,  gain  occurs  only  where 
the  is  dissociated  and  where  the  O^(^A)  fraction  is  greater  than 
0.15.  Unfortunately,  in  this  case,  much  of  the  I  atom  population  is  in 
the  O^(^A)  poor  jet  core.  This  illustrates  an  important  design 
consideration  for  COIL  devices.  Although  nonideal  mixing  may  result  in 
faster  dissociation  (c.f.  previous  section),  the  resulting  I  atoms  must 
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Figure  20.  Transverse  O  Mixing,  I  Dissociation,  and  Gain  Profiles  for  the 


be  well  mixed  with  the  0^(V)  stream  to  obtain  good  power  extraction. 

Note  how  the  gain  region  moves  outward  toward  the  0^  freestream  as 
the  jet  expands,  even  though  most  of  the  I  atom  concentration  is  in  the 
jet  core.  It  now  seems  clear  that  the  drop  off  in  the  streamwise 
direction  of  the  dissociation  rate  for  the  jet  cases  in  Figure  14  occurs 
because  the  in  the  jet  core  becomes  isolated  from  the  0^(V)  in  the 
freestream  by  a  layer  of  completely  dissociated  iodine  atoms.  is 
removed  rapidly  by  the  I  atom  ptjunplng  reaction,  Eq  (2.1),  as  it  diffuses 
toward  the  jet.  The  dissociation  rate  becomes  limited  by  the  rate  of 
transport  of  into  the  O^(^A)  stream. 

These  observations  suggest  that  the  most  efficient  laser  nozzle 
design  should  result  in  regions  of  high  concentrations  initially  to 
promote  dissociation  and  then  rapidly  mix  the  0^  and  I  streams  to  obtain 
good  extraction  characteristics.  This,  in  fact,  seems  to  be  the  reason 
why  the  the  transverse,  transonic  injection  nozzle  configuration 
(Figures  3  and  6a,  Chapter  2)  employed  on  the  most  recent  COIL  devices 
was  much  more  successful  than  other  configurations. 

We  now  compare  the  structure  of  the  3/1  and  1/1  jet-mixed  cases 
with  the  6/1  case  in  order  to  examine  the  effect  the  varying  velocity 
ratio  has  on  the  mixing  and  dissociation  processes.  In  order  to  do  so, 
we  need  to  first  establish  a  common  basis  for  comparison.  The  most 
obvious  would  be  to  compare  the  profiles  at  the  same  streamwise 
locations,  but,  as  in  comparisons  with  the  premixed  cases  discussed  in 
the  proceeding  section,  the  results  are  confused  by  the  differing  time 
scales  for  the  jets.  For  the  1/1  jet  case,  the  jet  flow  propagates 
approximately  a  factor  of  six  slower  than  the  flow  for  the  6/1  case. 
The  iodine  will  be  exposed  to  the  O^('a)  for  six  times  as  long  (time  and 
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distance  being  approximately  linearly  relatable,  c.f.  previous  section). 
For  any  chemical  rate  process,  time  is  clearly  an  important  variable. 
Based  on  this  argument,  the  transverse  structure  of  the  jets  are 
compared  at  similar  times  where  the  jet  velocity  is  employed  to  relate 
time  and  distance.  Figures  21,  22,  and  23  compare  the  transverse 

profiles  for  the  6/1  jet  at  x  -  6  cm  to  the  3/1  jet  at  x  -  3  cm  the  1/1 
jet  at  X  -  1  cm.  Admittedly,  scaling  the  distances  by  the  jet  velocity 
is  an  over  simplification  of  the  development  of  the  jet,  however,  this 
first  order  comparison  does  provide  some  insight. 

The  difference  in  the  total  0^  mass  fraction,  Eq  (5.12),  and  the 
O^(^A)  fraction,  Eq  (5.11),  profiles  for  the  lower  velocity  ratio  cases 
can  be  seen  in  Figure  21.  In  the  6/1  jet  case,  molecular  diffusion  is 
aided  by  convection  (entrainment)  due  to  the  large  velocity  difference 
producing  the  largest  0^  mass  fraction  on  the  jet  centerline.  Transport 
in  the  1/1  case  is  due  solely  to  diffusion  which  results  in  a  much  lower 
0^  concentration  in  the  jet  core.  For  all  of  these  cases,  very  little 
of  the  transported  toward  the  centerline  survives  the  combination 

of  the  dissociation  and  pumping  processes.  The  corresponding 

comparison  of  the  total  mass  fraction,  Eq  (5.10),  and  dissociation 
fraction,  Eq  (5.8)  is  shown  in  Figure  22.  Note  that  the  higher  initial 
mass  fraction  on  the  centerline  for  the  1/1  and  3/1  case  are  due  to  a 
lower  secondary  helium  flow  rate  for  the  slower  jets  (c.f.  Table  IV). 
The  exoansion  of  the  jets  for  all  three  cases  is  nearly  the  same. 
Nearly  all  of  the  which  diffuses  into  the  0^  primary  freestream  is 
dissociatt  d.  The  higher  entrainment  rates  for  the  6/1  and  3/1  jet 
velocity  ratios,  however,  lead  to  much  higher  dissociation  levels  in 
the  jet  core.  As  a  result,  the  total  dissociation  rate  is  faster  for 
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f  the  Transverse  0^  Mixing  Profiles  between  the  1/1 
X  =«  Icm,  the  3/1  Jet  Case  at  x  =  3cm,  and  the  6/1  Jet 


Figure  23.  Comparison  of  the  Transverse  Gain  Profiles  between  the  1/1  Jet 

Case  at  X  =  1cm,  the  3/1  Jet  Case  at  x  -  3cm,  and  the  6/1  Jet  Case 


the  higher  jet  velocity  ratio  cases  (c.f.  Figure  14  for  1/1  at  1  cm,  3/1 
at  3  cm,  and  6/1  at  6cm).  The  gain  profiles  for  these  three  cases, 
Figure  23,  show  an  interesting  trend.  The  gain  for  the  6/1  case  is  low 
even  though  the  total  dissociation  is  high.  This  is  true  because  more 
of  the  O^(^A)  near  the  mixing  region  was  consumed  in  creating  a  higher 
level  of  dissociation.  This  reinforces  the  point  that  the  must  be 
both  dissociated  and  thoroughly  mixed  to  produce  good  power  extraction 
qualities . 

Effect  of  Water  Quenching 

One  of  the  simplifying  assumptions  made  for  the  mixing  cases 
reviewed  in  the  preceding  section  was  that  no  H^O  vapor  was  included  in 
the  initial  conditions.  As  discussed  earlier,  this  eliminates  all  of 
the  quenching  (loss)  reactions  for  the  reduced  eleven  reaction  sec 
(Table  I,  Chapter  III)  employed  for  these  calculations.  In  order  to 
demonstrate  the  role  that  H^O  vapor  (higher  loss  rates)  has  on  the 
dissociation  process,  2%  H^O  vapor  (%  of  the  0^  flow)  was  added  to  the 
initial  conditions  for  the  6/1  case  and  the  associated  premixed  case. 
For  the  6/1  jet-mixed  case,  the  H^O  vapor  was  added  to  the  0^  freestream 
only.  This  is  physically  realistic  because  H^O  vapor  is  introduced  by 
the  BHP/Cl^  generator  in  actual  devices.  The  integrated  effect 
of  H^O  on  the  I^  dissociation  fraction,  Eq  (5  1),  for  both  of  these 
cases  is  seen  in  Figure  24.  As  expected,  H^O  interferes  with  the 
dissociation  process  (by  quenching  O^(^S)  and  I  ,  the  initiators,  and 
I^^ .  the  intermediate)  resulting  in  slower  dissociation  rates  for  both 
the  premixed  and  the  jet-mixed  cases.  The  H^O  appears  to  have  a  more 
detrimental  effect  in  the  jet  case.  This  is  probably  due  to  the  higher 
concentration  of  H^O  in  the  primary  freestream  than  in  the  premixed  case 
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remixed  and  Jet-Mixed  Cases. 


(i.e.,  H^O  is  spread  over  a  larger  volume  in  the  premixed  case). 

Elimination  of  H  0  vapor  from  the  0  flow  has  always  been  recognized  as 
2  2 

important  to  COIL  performance;  however,  this  was  attributed  to  the  fast 

★ 

quenching  rate  for  I  ,  the  lasing  species.  Figure  24  shows  that  H^O 
plays  an  important  role  in  retarding  the  I^  dissociation  process  as 
well . 

The  dissociation  efficiency  parameter,  Eq  (5.7),  is  plotted  in 
Figure  25  for  these  H^O  cases.  As  discussed  in  an  earlier  section,  the 
higher  loss  rate  due  to  the  H^O  quenching  reactions  results  in  less 
efficient  (higher  N)  dissociation.  Note,  the  impact  on  the  jet-mixed 
case  is  more  severe;  however,  this  is  due  to  the  higher  H^O 
concentration  in  the  primary  freestream  in  the  mixing  case.  With  H^O, 
the  value  of  the  efficiency  factor  begins  to  approach  the  often  used 
range  of  from  5  to  7  at  100%  dissociation.  Figures  17  and  25  raise  some 
important  questions  about  the  use  of  the  empirical  factor  N  for  COIL 
performance  prediction.  These  numerical  cases  point  out  that  N  is 
clearly  not  a  constant.  The  apparent  value  of  N  varies  with  both  the 
loss  rates  (^2*^  concentration)  and  the  dwell  time  between  iodine 
injection  and  power  extraction.  N  is,  therefore,  dependent  on  both 
device  geometry  and  operating  conditions.  The  practice  of  using  a  fixed 
N  (which  is  derived  from  a  demonstrated  device  design)  to  make  large 
design  excursions  is,  therefore,  extremely  suspect. 
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for  the  6/1  Premixed  and  Jet-Mixed  Cases. 


VI ■  Conclusion 


The  observations  made  from  the  numerical  results  of  the  previous 
section  provide  some  insight  into  the  nature  of  the  iodine  mixing  and 
dissociation  process  and  its  impact  on  the  performance  of  COIL  systems. 
One  common  belief  is  that  better  mixing  of  the  oxygen  and  iodine  streams 
will  result  in  faster,  more  efficient  dissociation.  This  is  true  only 
indirectly.  Comparison  of  the  one - dimens ional ,  premixed  cases  with  the 
jet -mixed  cases  in  Figure  14  clearly  showed  that,  due  to  hand  off  in 

rate  control  from  0^(  S)  and  0^(  A)  to  the  faster  I  - ►  I^*  mechanism, 

the  less  than  perfectly  mixed  jet  cases  resulted  in  higher  dissociation 
rates,  at  least  initially.  Conversely,  the  gain  profiles  in  Figures 
18-20  and  23  illustrate  the  need  to  achieve  complete  mixing  in  order  to 
obtain  high  gains  (and  power  extraction).  Obviously,  a  compromise 
design  is  required  which  results  in  regions  of  high  initial  iodine 
concentrations  to  initiate  the  dissociation  and  then  in  thorough  mixing 
of  oxygen  and  atomic  iodine  streams.  This  in  fact,  seems  to  be  the  case 
for  the  highly  successful  transverse  injection  nozzle  currently  being 
employed. 

As  should  be  expected,  the  magnitude  of  the  loss  mechanisms  (due  to 

H^O  in  this  case)  have  the  strongest  influence  on  the  efficiency  of  the 

dissociation  process.  Dissociation  is  slowed  by  the  quenching  of 
1  *  1 

0^(  E),  I  ,  and  1^'.  Mixing  plays  a  role  only  to  the  extent  that  more 

losses  are  incurred  as  the  rate  of  dissociation  is  limited  by  the  mixing 

of  the  two  streams.  Figures  17  and  25  illustrate  the  fact  that  the 

efficiency  of  the  combined  iodine  mixing  and  dissociation  process  as 

defined  by  N,  the  number  of  0  (^A)'s  required  to  dissociate  one  I  ,  is 

2  2 
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very  much  dependent  on  the  individual  device  design  and  operating 
conditions.  Larger  dwell  times  (the  flow  time  between  the  injection 
and  power  extraction)  and  greater  losses  (due  to  H^O)  result  in  a  higher 
apparent  N. 

The  recent  Roto-COIL  experience  discussed  in  Chapter  I  reinforces 
the  unacceptability  of  using  empirical  factors  like  N  to  predict  the 
performance  of  off-baseline  designs.  Experience  with  one - dimens ional , 
premixed  models  has  proven  to  be  just  as  unsatisfying.  As  demonstrated 
here,  the  premixed  results  display  a  very  different  character  altogether 
from  the  mixing  cases .  This  effect  may  become  even  more  pronounced  when 
attempting  to  model  the  performance  of  a  highly  three-dimensional  laser 
nozzle.  Clearly,  COIL  performance  prediction  will  require  a  more 
sophisticated  numerical  approach.  The  modeler,  however,  must  never 
forget  that  the  accuracy  of  any  prediction  is  only  as  valid  as  the  data 
base  supporting  it. 

The  most  basic  question  to  be  addressed  by  the  proposed 
experimental  and  computational  effort  is;  are  the  current  COIL  kinetics 
mechanism  and  rates  accurate?  Because  dissociation  of  iodine  in  an 
O^(^A)  stream  is  so  highly  dependent  on  the  mixing  as  well  as  the 
kinetics,  neither  the  simple  experiments  from  which  the  kinetics  data 
were  derived  nor  the  laser  device  experiments  have  been  able  to  provide 
an  answer  to  this  question.  The  planned  parallel  jet  experiment 
described  in  Appendix  A  is  intended  to  shed  some  light  on  this  subject; 
however,  just  as  the  results  of  an  unanchored  code  are  of  little  value, 
the  experimental  results  will  be  of  little  value  if  they  can  not  be 
accurately  modeled  using  the  code  developed  here  or  a  more  sophisticated 
one.  Several  issues  concerning  this  question  were  raised  lu  Chapter  V. 
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Certainly,  the  most  critical  is  whether  the  flow  will  remain  truly 
laminar.  The  modeling  difficulty  will  obviously  increase  dramatically 
if  this  is  not  so. 

The  nvimerical  experience  of  Chapter  V  does  provide  some  useful 
guidance.  One  effect  not  anticipated  in  the  design  of  the  experiment 
can  be  seen  in  Figure  14.  Note  that  although  the  comparisons  made  of 
the  dissociation  rates  for  the  three  jet  cases  based  on  time  did  provide 
some  insight,  th“  spatial  profiles  display  very  similar  character.  This 
will  make  diagnosing  the  experimental  flows  much  more  difficult.  On  a 
more  promising  note,  the  gain  profiles  in  Figures  18-20  and  23  show  a 
much  more  distinct  behavior.  Gain  mapping  is  a  promising  diagnostic 
technique  because  it  provides  information  not  only  on  the  dissociation, 
but,  also  on  the  mixing  of  the  atomic  iodine  with  the  0^(V). 

The  most  immediate  goals  of  the  proposed  experimental  and 
computational  effort  are  to  acquire  a  strong  dissociation  data  base  and 
validate  the  current  AFWL  COIL  kinetics  data  set.  The  ultimate  success 
of  the  combined  experimental  and  computational  effort  lies  in 
identifying  new  techniques  for  determining  fundamental  rate  information. 
This  combined  approach  is  not  so  different  from  the  techniques  employed 
to  measure  the  original  rate  data  (9)  (11).  Analytical  models  for 
possible  mechanisms  were  compared  to  overall  creation  and  depletion 
rates  measured  on  simple  flow  tube  experiments.  Figure  5,  Chapter  II, 
points  to  some  possible  difficulties  with  this  method  which  may  have 
resulted  from  attempting  to  oversimplify  the  problem.  It  remains  to  be 
shown  that  all  of  the  Important  influences  can  indeed  be  modeled  in  the 
more  complicated  experiments. 
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Appendix  A 

Concurrent  Experiment:  Configuration  and  Expected  Operating  Conditions 

Figure  Al  is  a  schematic  of  the  COIL  device  which  will  be  used  in 
planned  experiments.  The  major  components  of  the  excited  oxygen 
delivery  system  include  a  BHP/CI2  sparger,  liquid  separator,  transition 
and  vacuum  bypass,  H^O  cold  trap.  The  nozzle  cavity  allows  optical 
diagnostic  access  to  both  the  side  and  top  of  the  two  dimensional  slit 
injector.  Pressure  taps  spaced  at  half  inch  intervals  along  the  top  and 
bottom  walls  allow  for  an  experimental  measurement  of  the  streamwise 
pressure  gradient.  The  initial  conditions  of  the  oxygen  stream  are 
measured  in  a  diagnostic  duct  located  just  upstream  of  the  slit 
injector.  Iodine  is  supplied  by  passing  heated  helium  over  a  heated 
iodine  bed.  The  concentration  of  the  entrained  iodine  is  measured 
optically  just  before  Injection.  The  transverse  injection  supersonic 
slit  nozzle  shown  in  Figure  6a  exists  and  can  be  mated  to  the  same 
O^(^A)  delivery  system  for  comparison  with  the  two  dimensional  results. 

Figure  A2  shows  the  detailed  geometry  of  the  two-dimensional  iodine 
nozzle  and  the  optical  cavity.  The  cross  section  of  the  test  section  is 
1.25  X  2.0  inches  (3.175  X  5.08  cm)  The  height  of  the  two  dimensional 
nozzle  is  0.075  Inches  (0.190  cm).  The  baseline  flow  rates  through  the 
bH?/Cl^  reactor  are  30  mmole/sec  of  chlorine  at  a  3/1  diluent  ratio  (90 
mmole/sec  of  helium)  and  100%  chlorine  utili^;atior.  (complete  conversion 
of  chlorine  to  oxygen).  The  generator  is  expected  to  deliver  an  excited 
oxygen  fraction  (0^(‘A)/0^(total) )  of  0.40.  The  temperature  of  the 
oxygen  stream  will  be  approximately  300  K  as  it  exits  the  cold  trap. 
The  nominal  iodine  to  oxygen  ratio  is  0.02  (0.6  mmole/sec  of  I^) .  The 
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Diagnostic  Duct  Diagnos 
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Figure  Al.  COIL  Diagnostic  Flow  Tube  Assembly 
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Figure  A2 .  COIL  Diagnostic  Flow  Tube  Injection  Nozzle 


velocity  of  the  iodine  jet  is  controlled  by  varying  the  secondary 
diluent  (helium  in  this  case).  The  nozzle  was  designed  for  a  6/1 
velocity  ratio  for  a  secondary  diluent  flow  rate  of  30  mmole/sec  of 
helium.  The  iodine  flow  is  heated  to  400  K  to  avoid  condensing  iodine 
vapor  on  the  flow  surfaces.  The  pressure  in  the  cavity  is  set  by  a 
0.25  X  2.0  inch  (0.635  X  5.08  cm)  choke  plate  located  well  downstream  of 
the  nozzle  exit  plane.  The  cavity  pressure  is  approximately  25  torr 
(3.33x10  dyne/cm).  The  velocity  of  the  freestream  at  these  flow 
conditions  is  5,550  cm/sec  (Ms  0.1)  and  the  velocity  of  the  jet  is 
33,000  cm/sec  (M  s  0.4). 

The  basic  array  of  diagnostics  for  this  experiment  include  a  number 
of  flow  and  media  diagnostics  of  increasing  complexity.  Flow 
visualization  using  one  of  two  techniques  will  be  used  to  provide  a 
qualitative  picture  of  the  flow  field  and  mixing  region.  The  first 
technique.  Laser  Induced  Fluorescence  (LIF) ,  uses  an  Argon  Ion  (Ar^) 
laser  to  excite  an  electronic  transition  in  the  I2  molecule  during  cold 
flows  (no  02(^4)). 

I2(X)  +  5145nm  - >  12(B)  (A.l) 

This  same  transition  is  also  pumped  in  the  02(^A)/l2  system  giving  off  a 

characteristic  yellow  flame.  Infrared  spectroscopy  using  spatially 

filtered  intrinsic  germanium  photo  diode  detectors  will  be  used  to 
1  ^ 

monitor  02(  A)  and  I  species  concentrations.  Residual  (unreacted) 
chlorine  is  measured  by  absorption  using  a  Helium-Cadmium  laser  in  the 

diagnostic  duct.  I  and  I^  can  also  be  measured  by  absorption  using  an 

*  + 

I  probe  laser  (or  diode  laser)  and  an  Ar  laser,  respectively. 
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Appendix  B 

Coefficient?  for  Generalized  Finite-Difference  Operators 

2  2 

The  general  finite -difference  operators  for  df/d-ip  and  3  f/dip 
defined  in  Chapter  IV,  Eqs  (4.8)  and  (4.9)  are  repeated  here  for 
convenience . 


f(x.^) 


D 


(B.l) 


^2  f(x.^) 


A  f(x.^  )  +  B  f(x.^  )  +  C  f(x,^  ) 

Z  A  Z  Z  U 


(B.2) 


The  values  for  the  coefficients  A,  B,  C,  A,  B,  C,  D,  i>  ,  ip  ,  and 

1  i’  i’  2’  2’  2  ’ 

iP^  are  given  in  Table  B.  The  notation  used  is  shown  to  the  right  of  the 
schematics  in  Figure  7,  Chapter  IV. 

TABLE  B 

Coefficients  for  the  Generalize  Finite-Difference  Operators 
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1 
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PROGRAM  SPECIES  WRITTEN  BY; 

JOHN  S.  LAWRENCE,  GAE-83D 

MIXING  WRITTEN  BY: 

CAPT.  DAVID  BJURSTROM 
CAPT.  DENISE  PEREZ 
LT.  JEFFREY  MILLER 
GAE  84D 

COILFT  WRITTEN  BY: 

CAPT  JEFFREY  A.  MILLER 
MASTERS  MARCH,  1989 


THE  PURPOSE  OF  THIS  PROGRAM  IS  TO  CALCULATE  THE  PROPERTIES  OF  A 
MULTI -SPECIES  FLOW  AS  IT  PROPAGATES  DOWNSTREAM,  TAKING  INTO  ACCOUNT 
TRANSPORT  PHENOMENA  DUE  TO  MOMENTUM,  ENERGY,  AND  MASS  OF  SPECIES 
(INCLUDING  CHEMICAL  REACTION).  THE  USER  MUST  SPECIFY  THE  INITIAL 
VELOCITY  AND  TEMPERATURE  PROFILES  AND  THE  INITIAL  MOLE  FRACTIONS 
OF  EACH  SPECIES. 

IN  USER'S  REFERENCES  THROUGHOUT  THE  PROGRAM,  THE  FOLLOWING 
ABBREVIATIONS  APPLY: 

BS&L  -  BIRD,  STEWART,  AND  LIGHTFOOT,  TRANSPORT  PHENOMENA,  1960. 
V&K  -  VINCENTI  AND  KRUGER,  PHYSICAL  GAS  DYNAMICS,  1965. 

JANZ  -  JANZ,  THERMODYNAMIC  PROPERTIES  OF  ORGANIC  COMPOUNDS,  1967. 

SYMBOLOGY  WITHIN  THE  PROGRAM  IS  AS  FOLLOWS: 

X  -  POSITION  VARIABLE  TANGENT  TO  THE  FREE-STREAM  FLOW.  (CM) 

DX  -  INCREMENT  IN  THE  X-DIRECTION. (CM) 

Y  -  POSITION  VARIABLE  NORMAL  TO  THE  FREE-STREAM  FLOW.  (CM) 

DY  -  INCREMENT  IN  THE  Y- DIRECTION.  (CM) 

PSI  -  POSITION  VARIABLE  NORMAL  TO  THE  STREAMWISE  VELOCITY  (PSI). 
(G/CM*S) 

K  -  INCREMENT  BETWEEN  LINES  OF  CONSTANT  PSI.  (G/CM*S) 

I  AND  ICOUNT  -  INTEGER  COUNTERS  IN  THE  X-DIRECTION. 

J  AND  L  -  INTEGER  COUNTERS  IN  THE  Y- DIRECTION. 

M  AND  N  -  NUMBER  OF  DESIRED  STEPS  IN  THE  X-  AND  Y-DIRECTIONS . 

NS  -  NUMBER  OF  SPECIES. 

Q  AND  S  -  INTEGER  COUNTERS  FOR  SPECIES. 

U  -  VELOCITY.  (CM/S) 

DIU  AND  D2U  -  1ST  AND  2ND  DERIVATIVES  OF  U  WRT  PSI. 

T  -  TEMPERATURE.  (K) 

DIT  AND  D2T  -  1ST  AND  2ND  DERIVATIVES  OF  T  WRT  PSI. 
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*  PR  -  PRESSURE.  (DYNE/CM**2) 

*  DPDX  -  PRESSURE  DERIVATIVE  WRT  X. 

*  R  -  DENSITY  (RHO).  (G/CM**3) 

*  MU  -  VISCOSITY.  (G/CM*S) 

*  MUS  -  VISCOSITY  OF  SPECIES.  (G/CM*S) 

*  KT  -  THERMAL  CONDUCTIVITY.  (DYNE/S*K) 

*  KTS  -  THERMAL  CONDUCTIVITY  OF  SPECIES.  (DYNE/S*K) 

*  CP  -  SPECIFIC  HEAT  AT  CONSTANT  PRESSURE.  (DYNE*CM/K) 

*  CPS  -  SPECIFIC  HEAT  OF  SPECIES.  (DYNE*CM/MOLE*K) 

*  H  AND  HT  -  ENTHALPY  OF  SPECIES.  (DYNE*CM/G) 

*  DIH  AND  DIHT  -  1ST  DERIVATIVE  OF  ENTHALPY  WRT  PSI. 

*  DIF  -  DIFFUSIVITY  OF  SPECIES.  (CM**2/S) 

*  DIFS  -  BINARY  DIFFUSION  COEFFICIENT.  (CM**2/S) 

*  MW  -  MOLECULAR  WEIGHT.  (GM/MOLE) 

*  RBAR  AND  RCA  -  UNIVERSAL  GAS  CONSTANT.  (DYNE*AM/MOLE*K)  OR 

*  (CAL/MOLE*K) 

*  EPSK  AND  SIG  -  LENNARD-JONES  PARAMETERS  (EPSILON/K  AND  SIGMA). 

*  OM  -  CHAPMAN- ENSKOG  PARAMETER  (OMEGA). 

*  MRU  AND  KRU  -  MUXRXU  AND  KTXRXU,  TREATED  AS  SEPAR.\TE  VARIABLES. 

*  DIMRU  AND  DIKRU  -  CORRESPONDING  1ST  DERIVATIVES  WRT  PSI. 

*  W  AND  WT  -  MASS  FRACTION  OF  SPECIES. 

*  DIW,  D2W,  DIWT,  AND  D2WT  -  1ST  AND  2ND  DERIVATIVES  OF  W  WRT  PSI. 

*  MF  -  MOLE  FRACTION  OF  SPECIES. 

*  WD  -  CREATION  OF  MASS  (RATE)  TERM  FROM  CHEMICAL  REACTION. 

*  MWM  -  MOLECULAR  WEIGHT  OF  THE  MIXTURE.  (G/MOLE) 

*  PHI  -  SPECIES  TRANSPORT  COEFFICIENT. 

*  XPHI  -  SUMMATION  OF  MFXPHI  OVER  SPECIES. 

*  GRU  AND  GRT  -  DIFX(R**2)XU,  TREATED  AS  A  SEPARATE  VARIABLE. 

*  DIGRU  -  CORRESPONDING  1ST  DERIVATIVE  WRT  PSI. 

*  TTABLE  -  TEMPERATURE  TABLE  FOR  CP  AND  H  INTERPOLATION.  (K) 

*  CPTBLE  -  CP  TABLE  ENTRIES.  (CAL/MOLE*K) 

*  HTABLE  -  H  TABLE  ENTRIES.  (KCAL/MOLE) 

*  C  -  CONSTANT  0.0<C<-1.0.  USED  FOR  DETERMINATION  OF  DX. 

* 

INTEGER  I,J,M,N,ICOUNT,S,NS 
PARAMETER  (NS-10,N-32) 

REAL*4  C 

C  PARAMETER  (C-0.05) 

R£AL*4  RBAR 

PARAMETER  (RBAR-8 . 3143E+07) 

REAL*8  U(0:N) ,R(0:N) ,T(0:N) ,W(NS.0:N) 

REAL*8  UIP1(0;N) ,TIP1(0:N) ,WIP1(NS.0:N) 

REAL*8  WD(NS,0:N) 

REAL*8  X,DX,DXP,PR,DPDX,TIME(0:N) 

REAL*8  Y(0;N),DY(N),K(N),PSI(0:N) 

REAL*8  MW.MWM.MF 
REAL*8  SIG, EPSK 
REAL*8  A.B.HFORM 

REAL*8  MU(0;N) ,CP(0;N) ,H(NS,0:N) ,KT(0:N) ,DIF(NS,0;N) 

INTEGER  IA(0:N) ,IB(0:N) ,IC(0:N) 

REAL*8  A1(0;N) ,B1(0:N) ,C1(0:N) 

REAL*8  A2(0:N) ,B2(0:N) ,C2(0:N) ,D(0:N) 

REAL*8  D1U,D2U, MRU, DIMRU 
REAL*8  D1T,D2T, KRU, DIKRU 
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REAL*8  D1W.D2W,D1H,GRU.D1GRU 
REAL*8  TERM1,TERM2,DEN0M 
INTEGER  SSAVE.JSAVE 

COMMON  /PROPl/  MW(NS) ,MWM(0:N) ,MF(NS,0;N) 

COMMON  /LENJONES/  SIG(NS) , EPSK(NS) 

COMMON  /SPHTl/  A(NS , 0 ; 4) , B(NS ,0 :4) ,HFORM(NS) 

* 

*  ESTABLISH  INPUT  (MIXIN)  AND  OUTPUT  (MIXOUT)  FILES.  THE  INPUT  FILE 

*  IS  UNFORMATTED  ,  AND  SHOULD  BE  CONSTRUCTED  AS  FOLLOWS: 

*  FIRST  LINE 

*  M  -  NUMBER  OF  STREAMWISE  STEPS 

*  C  -  STABILITY  FACTOR 

*  SECOND  LINE 

*  X  -  INITI vL  X 

*  PR  -  PRESSURE 

*  DPDX  -  PRESSURE  GRADIENT 

*  NEXT  NS*3  LINES 

*  (MW,  ESPK,  SIG)  -  MOLECULAR  WEIGHT.  LENNARD- JONES  POTENTIALS 

*  (A(S.I) ,1-0,4)  -  COEFFICIENTS  FOR  FOURTH  ORDER  POLYNOMIAL  FITS  FOR 

*  (B(S,I) ,1-0,4)  -  CP  AND  H.  DATA  FROM  JANAF  TABLES 

*  NEXT  N  LINES  -  ONE  LINE  ENTRY  FOR  EACH  PS I  STEP  STARTING  AT  THE 

*  BOTTOM  BOUNDRY.  EACH  LINE  CONTAINS  Y(I) ,U(I) ,T(I) ,MF(1 , I) - 

*  MF(NS,I). 

* 

OPEN(UNIT-15 , FILE-' coilin. dat ' ) 

REWIND  15 

OPEN (UNIT-16 , FILE- ' coilout . dat ' ) 

REWIND  16 

* 

NMl-N-1 
PSI(0)-0.0 
SSAVE-O 
J SAVE-0 

* 

*  EXTRACT  DATA  FROM  INPUT  FILE. 

* 

READ(15,*)  M,C 
READ(15,*)  X, PR, DPDX 
READ(15,*)  IPRINT.XPRINT.DXPRNT 
DO  5  S-1,NS 

READ(15,*)  MW(S) ,EPSK(S),SIG(S) 

READ(15,*)  (A(S,I),I-0.4) 

READ(15,*)  (B(S,I) ,1-0,4) ,HFORM(S) 

5  CONTINUE 

DO  10  J-0,N 

READ(15,*)  Y(J),U(J),T(J).(MF(II.J),II-1,NS) 

Y(J)-Y(J)*2.54 
10  CONTINUE 

* 

*  DENSITY  AND  MASS  FRACTIONS  ARE  COMPUTED  BY  FIRST  OBTAINING  THE 

*  MOLECULAR  WEIGHT  OF  THE  MIXTURE  (MWM  -  SUMMATION  OF  MFXMW  OVER 

*  ALL  SPECIES),  COMPUTING  DENSITY  OF  THE  MIXTURE  BY  THE  PERFECT  GAS 

*  LAW,  AND  TAKING  THE  APPROPRIATE  RATIO  FOR  MASS  FRACTION.  SEE 

*  BS&L,  TABLE  16.1-1,  P498. 
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DO  25  J-0,N 
MWM(J)-0.0 
DO  15  S-1,NS 

MWM ( J ) -MWM ( J ) +M F ( S , J ) *MW ( S ) 

15  CONTINUE 

R( J ) -MWM( J ) *PR/ (RBAR*T ( J ) ) 

DO  20  S-1,NS 

W(S,J)-MF(S,J)*MW(S)/MWM(J) 

20  CONTINUE 

25  CONTINUE 

* 

*  THE  VON  MISES  TRANSFORMATION,  D(PSI)/Dy-RHO*U,  IS  USED  TO  COMPUTE 

*  VALUES  FOR  PS I  AND  K,  WHERE  RHO  AND  U  ARE  AVERAGED  OVER  EACH  STEP 

*  OUTWARD  FROM  THE  WALL,  AND  THE  DERIVATIVE  IS  APPROXIMATED  BY  THE 

*  RATIO  OF  INCREMENTAL  VALUES. 

* 

DO  30  J-1,N 

DY(J)-Y(J)-Y(J-1) 

K(J)-DY(J)*(R(J)+R(J-l))*(U(J)+U(J-l))/4. 

PS1(J)-PSI(J-1)+K(J) 

30  CONTINUE 

* 

*  SET  GENERALIZED  FINITE-DIFFERENCE  COEFFICIENTS.  WALL  BOUNDARY 

*  CONDITIONS  ARE  COMMENTED  OUT  IN  FAVOR  OF  REFLECTION  BOUNDARY 

*  CONDITIONS . 

* 

DO  31  J-0,N 
IF(J.EQ.O)THEN 
lA(0)-2 
IB(0)-1 
IC(0)-0 

C  A1(0)— K(1)*K(1) 

C  B1(0)-(K(1)+K(2))*(K(1)+K(2)) 

C  Cl(0)-(-Al(0)-Bl(0)) 

C  A2(0)-K(l) 

C  B2(0)— (K(l)+K(2)) 

C  C2(0)-K(2) 

C  D(0)-K(1)*K(1)*K(2)+K(1)*K(2)*K(2) 

Al(0)-0.0 

Rl(0)-0.0 

31(0)-0.0 

A2(0)-0.0 

B2(0)-1.0 

C2(0)— 1.0 

D(0)-K(1)*K(1) 

ELSE  IF(J.EQ.N)THEN 
lA(N)-N-2 
IB(N)-N-1 
1C(N)-N 
Al(N)-0.0 
Bl(N)-0.0 
Cl(N)-0.0 
A2(N)-0.0 
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B2(N)-1.0 
C2(N)  — 1  .0 
D(N)-K(N)*K(N) 

C  Al(N)-K(N)*K(N) 

C  B1(N)— (K(N)+K(N-1))*(K(N)+K(N-1)) 

C  C1(N)— (A1(J)+B1(J)) 

C  A2(N)-K(N) 

C  B2(N)— (K(N)+K(N-1)) 

C  C2(N)-K(N-1) 

C  D(N)-(K(N)*K(N)*K(N-1)+K(N)*K(N-1)*K(N-1) ) 

ELSE 

IA(J)-J+1 

IB(J)-J 

IC(J)-J-1 

A1(J)-K(J)*K(J) 

B1(J)— (K(J)*K(J) -K(J+1)*K(J+1) ) 

C1(J)  — K(J+1)*K(J+1) 

A2(J)-K(J) 

B2(J)— (K(J)+K(J+1)) 

C2(J)-K(J+1) 

D(J)-(K(J)*K(J)*K(J+1)+K(J)*K(J+1)*K(J+1)) 

END  IF 
31  CONTINUE 

★ 

*  COMPUTE  INITIAL  VISCOSITY,  SPECIFIC  HEAT,  ENTHALPY,  THERMAL 

*  CONDUCTIVITY,  DIFFUSIVITY,  AND  MASS  CREATION  TERM 

* 

CALL  PR0P(PR,T,MU,CP,H,KT,DIF) 

CALL  RATE(W,R.T.WD,MW) 

* 

*  RECORD  INITIAL  PROPERTIES  OF  THE  FLOW 

* 

I-O 

WRITE(16,200)  N,NS,M,C 
WRITE(16,210)  I.X.DX.PR.DPDX 
DO  35  J-0,N 
TIME(J)-0.0 

WRITE(16,220)  PSI(J) ,Y(J) .TIME(J) ,R(J) ,U(J) ,T(J) .MWM(J) , 

+  MU(J),CP(J).KT(J) 

35  CONTINUE 

DO  36  J-0,N 

VRITE(16,220)  (MF(S , J) , S-1 ,NS) 

36  CONTINUE 

DO  37  J-0,N 

WRITE(16,220)  (WD(S , J) , S-1 . NS) 

37  CONTINUE 

DO  38  J-0,N 

WRITE(16,220)  (H(S , J) . S-1 ,NS) 

38  CONTINUE 

no  39  J-0,N 

WRITE(16,220)  (DIF(S , J) . S-1 ,NS) 

39  CONTINUE 

* 

*  THIS  LOOP  STEPS  THE  ENTIRE  SOLUTION  FORWARD  BY  AN  INCREMENT  IN  THE 
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*  X- DIRECTION. 

* 

lCOUNT-0 
DO  85  I-1,M 

* 

*  ESTABLISH  DESIRED  SIZE  OF  INCREMENTAL  STEP  IN  X-DIRECTION. 

* 

DX-0 . 1 

C  IF(X.GT.1.0)C-0.05 

ICOUNT-ICOUNT+1 
DO  45  S-1,NS 
DO  40  J-1,NM1 

AMPH-(DIF(S,J)*R(J)**2*U(J)+DIF(S,J+1)*R(J+1)**2 
+  *U(J+l))/2.0 

AMMH-(D1F(S,J)*R(J)**2*U(J)+DIF(S,J-1)*R(J-1)**2 
+  *U(J-l))/2.0 

IF(AMPH.NE.O.O.AND.AMMH.NE.O.O)THEN 

DXP-C*(K(J)*K(J+1)*(K(J)+K(J+1)))/(2.0*(K(J)*AMPH 
+  +K(J+1)*AMMH)) 

IF(DXP.LT.DX)THEN 
SSAVE-S 
JSAVE-J 
DX-DXP 
END  IF 
END  IF 

40  CONTINUE 

45  CONTINUE 

DO  55  J-0,N 

* 

*  COMPUTE  NEW  VALUES  FOR  VELOCITY,  TEMPERATURE,  AND  MASS  FRACTION. 

*  THESE  DIFFERENCE  EQUATIONS  ARE  DERIVED  FROM  A  CONTROL  VOLUME 

*  ANALYSIS  OF  CONTINUITY,  MOMENTUM,  AND  ENERGY,  WITH  BOUNDARY  LAYER 

*  ASSUMPTIONS.  SEE  ACCOMPANYING  DERIVATIONS. 

* 

* 

*  SOLVE  THE  MOMENTUM  EQUATION. 

★ 

MRU-MU(J)*R(J)*U(J) 

D1U-(A1(J)*U(IA(J))+B1(J)*U(IB(J))+C1(J)*U(IC(J)))/D(J) 

D1MRU-(A1(J)*MU(IA(J))*R(IA(J))*U(IA(J))+ 

+  B1(J)*MU(IB(J))*R(IB(J))*U(IB(J))+ 

+  C1(J)*MU(IC(J))*R(IC(J))*U(IC(J)))/D(J) 

D2U-(A2(J)*U(IA(J))+B2(J)*U(IB(J))+C2(J)*U(IC(J)))*2./D(J) 
UIPl ( J ) -U ( J ) +DX* ( - DPDX/ (R( J ) *U ( J ) ) +D1MRU*D1U+MRU*D2U ) 

* 

*  SOLVE  EACH  OF  THE  SPECIES  EQUATIONS  AND  COMPUTE  TERMS  OF  THE  ENERGY 

*  EAUTION. 

* 

TERMl-0 . 0 
TERM2-0 . 0 
DO  50  S-1,NS 

GRU-DIF(S,J)*U(J)*R(J)**2 

D1H-(A1(J)*H(S,IA(J))+B1(J)*H(S,IB(J))+ 

+  C1(J)*H(S,IC(J)))/D(J) 
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D1W-(A1(J)*W(S,IA(J))+B1(J)*W(S,IB(J))+ 

+  C1(J)*W(S.IC(J)))/D(J) 

D1GRU-(A1(J)*D1F(S,1A(J))*R(IA(J))**2*U(IA(J))+ 

+  B1(J)*DIF(S,IB(J))*R(IB(J))**2*U(IB(J))+ 

+  C1(J)*DIF(S,IC(J))*R(IC(J))**2 

+  *U(IC(J)))/D(J) 

D2W-(A2(J)*W(S,1A(J))+B2(J)*W(S,IB(J))+ 

+  C2(J)*W(S,IC(J)))*2./D(J) 

WIP1(S,J)-W(S,J)+DX*(WD(S,J)/(R(J)*U(J))+D1GRU*D1W 
+  +GRU*D2W) 

TERMl-TERMl+WD ( S , J ) *H ( S , J ) 

TERM2-TERM2+GRU*D1W*D1H 
50  CONTINUE 

SOLVE  THE  ENERGY  EQUATION. 

KRU-KT(J)*R(J)*U(J) 

D1T-(A1(J)*T(IA(J))+B1(J)*T(IB(J))+C1(J)*T(IC(J)))/D(J) 

D1KRU-(A1(J)*KT(IA(J))*R(IA(J))*U(IA(J))+ 

+  B1(J)*KT(IB(J))*R(IB(J))*U(IB(J))+ 

+  C1(J)*KT(IC(J))*R(IC(J))*U(IC(J)))/D(J) 

D2T-(A2(J)*T(IA(J))+B2(J)*T(IB(J))+C2(J)*T(IC(J)))*2./D(J) 

TIP1(J)-T(J)+DX*(DPDX/(R(J)*CP(J))+MRU*D1U**2/CP(J) 

+  +(D1KRU*D1T+KRU*D2T)/CP(J)) 

TIP1(J)-TIP1(J) -DX*(TERM1/(R(J)*U(J)*CP(J) ) -TERM2/CP(J) ) 
TIME(J)-TIME(J)+DX/U(J) 

55  CONTINUE 

X-X+DX 

THE  MOLE  FRACTIONS  ARE  COMPUTED  BY  TAKING  THE  RATIO  OF  W/MW 
TO  THE  SUMMATION  OF  W/MW  OVER  ALL  SPECIES.  THE  DENSITY  IS 
THEN  COMPUTED  IN  THE  SAME  MANNER  AS  PREVIOUSLY,  AND  Y  AND  DY  ARE 
BACKED  OUT  OF  THE  SAME  EQUATION  USED  TO  CALCULATE  PS I  AND  K. 

DO  70  J-0,N 
U(J)-UIP1(J) 

T(J)-TIP1(J) 

DENOM-0.0 
DO  60  S-1,NS 

W(S,J)-WIP1(S,J) 

DENOM-DENOM+W ( S , J ) /MW ( S ) 

60  CONTINUE 

MWM(J)-l./DENOM 
DO  65  S-1,NS 

MF(S,J)-W(S,J)*MWM(J)/MW(S) 

65  CONTINUE 

R( J ) -MWM( J ) *PR/ (RBAR*T ( J ) ) 

70  CONTINUE 

DO  75  J-1,N 

DY(J)-4.*K(J)/((R(J)+R(J-1))*(U(J)+U(J-1))) 

Y(J)-Y(J-1)+DY(J) 

75  CONTINUE 

COMPUTE  VISCOSITY,  SPECIFIC  HEAT,  ENTHALPY,  THERMAL  CONDUCTIVITY, 


*  DIFFUSIVITY,  AND  MASS  CREATION  TERM. 

* 

CALL  PR0P(PR,T,MU,CP,H,KT,D1F) 

CALL  RATE(W,R,T,WD,MW) 

* 

*  RECORD  PROPERTIES  OF  THE  FLOW.  THE  ICOUNT  PROCEDURE  IS  DETERMINED 

*  BY  HOW  OFTEN  THE  USER  WISHES  A  PRINT  OUT  OF  RESULTS. 

* 

IF( ICOUNT . EQ . IPRINT . OR . X . GT . XPRINT ) THEN 
ICOUNT-0 

IF(X.LT.1.0)THEN 

XPRINT-XPRINT+DXPRNT/2 . 0 
ELSE 

XPRINT-XPRINT+DXPRNT 
END  IF 

WRITE(16,210)  I,X,DX,PR.DPDX 
DO  80  J-0,N 

WRITE(16,220)  PSI(J) ,Y(J) ,TIME(J) ,R(J) ,U(J) ,T(J) .MWM(J) , 

+  MU(J),CP(J).KT(J) 

80  CONTINUE 

DO  81  J-0,N 

WRITE(16,220)  (MF(S,J) .S-1,NS) 

81  CONTINUE 

C  DO  82  J-0,N 

C  WRITE(16,220)  (WD(S,J) ,S-1,NS) 

C  82  CONTINUE 

C  DO  83  J-0,N 

C  WRITE(16,220)  (H(S . J) ,S-1 ,NS) 

C  83  CONTINUE 

C  DO  84  J-0,N 

C  WRITE(16,220)  (DIF(S,J) ,S-1,NS) 

C  84  CONTINUE 

END  IF 
85  CONTINUE 

200  FORMATC  '  ,  3(110) ,  E13 . 5) 

210  FORMATS  M5,4(E13.5)) 

220  FORMATC  MO(E13.5)) 

STOP 

END 

*'**-**'******'*r**-.':***-*'******'<r*****************************************Sk’*** 


*  * 

*  THIS  SUBROUTINE  COMPUTES  THE  THERMODYNAMIC  PROPERTIES  AND  THE  * 

*  TRANSPROT  COEFFICIENTS.  * 

*  * 


************************************************************************** 

* 

*  NOTE:  1  (CAL)  -  4.186E+07  (DYNE*CM) 

*  1  (DYNE/CM**2)  -  9.869E-07  (ATM) 

*  RBAR  -  8.3143E+07  (ERG/MOLE*K) 

"k 

SUBROUTINE  PROP(PR,T,MU,CP,H,KT,DIF) 

INTEGER  L,N,q,NS,QS 
PARAMETER  (NS-10,N-32) 

REAL*4  RBAR 
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PARAMETER  (RBAR-8 . 3143E+07) 

REAL*8  T(0:N) , PR, PRESS 

REAL*8  MU(0:N),CP(0:N).H(NS.0:N),KT(0:N) ,DIF(NS,0;N) 

REAL*8  MW,MWM,MF 

REAL*8  SIG.EPSK 

REAL*8  A.B.HFORM 

REAL*8  MUS (NS ) . OMEGA , PHI . XPHI 

REAL*8  CPS.KTS 

REAL*8  D1FS,DIFSUM,EPSKT,SIGT 

COMMON  /PROPl/  MW(NS) ,MWM(0;N) ,MF(NS,0:N) 

COMMON  /LENJONES/  SIG(NS) ,EPSK(NS) 

COMMON  /SPHTl/  A(NS,0:4) ,B(NS,0;4) .HFORM(NS) 

* 

PRESS-PR*9.869E-07 
DO  20  L-0,N 
MU(L)-0.0 
CP(L)-0.0 
KT(L)-0.0 

* 

*  COMPUTE  THE  VISCOSITY  OF  THE  MIXTURE  AND  VISCOSITIES 

*  OF  THE  INDIVIDUAL  SPECIES  BY  EMPLOYING  THE  CHAPMAN -ENSKOG  FORMULAS. 

*  THE  VALUE  FOR  OMEGA  IS  DETERMINED  FROM  AN  ALGEBRAIC  CURVE- FIT 

*  ANALYSIS  OF  TABULATED  DATA  FOR  OMEGA  AS  A  FUNCTION  OF  EPSK  AND  T. 

*  SEE  BS&L;  EQNS  1,4-18,  1.4-19,  1.4-20,  PP22-25,  AND  TABLES  B1  AND 

*  B2,  PP744-746. 

* 

DO  5  Q-1,NS 

OMEGA-. 588-H./((l./EPSK(Q))*T(L)+. 125)**. 6666 
MUS(Q)-.000026693*SQRT(MW(Q)*T(L))/(OMEGA*SIG(Q)**2) 

5  CONTINUE 

DO  15  Q-1,NS 
XPHI-O . 0 
DIFSUM-O.O 
DO  10  QS-1,NS 

PHI-((1.+((MUS(Q)/MUS(QS))**.5)*(MW(QS)/MW(Q))**.25) 

+  **2)/(8.*(l.+MW(Q)/MW(QS)))**.5 

XPHI-XPHI+MF(QS ,L)*PHI 
IF(Q.EQ.QS)THEN 
DIFS-O.O 
ELSE 

EPSKT-SQRT(EPSK(Q)*EPSK(QS) ) 

SIGT-(SIG(Q)+SIG(QS))/2. 

OMEGA-. 488+1. /((l./EPSKT)*T(L)+. 125)**. 6666 
DIFS-.0018583*SQRT((1./MW(Q)+1./MW(QS))*T(L)**3)/ 

+  (PRESS*0MEGA*SIGT**2) 

DIFSUM-DIFSUM+MF(QS , L)/DIFS 
END  IF 

10  CONTINUE 

MU ( L) -MU ( L) +MF ( Q , L) *MUS (Q ) /XPHI 

* 

*  COMPUTE  THE  SPECIFIC  HEAT  AT  CONSTANT  PRESSURE  OF 

*  VALUES  FOR  THE  SPECIES  ARE  LINEAR  INTERPOLATED  BETWEEN  INPUT  VALUES 

*  SUPPLIED  BY  A  JANAF  TABLE.  NOTE:  UNITS  FOR  INPUT  ARE  CP  (CAL/MOLE* 

*  K)  AND  H(KCAL/MOLE) .  THESE  UNITS  ARE  CHANGED  TO  CGS  UNITS  IN  THE 
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*  SUBROUTINE.  CP  OF  THE  MIXTURE  IS  A  SUM  OF  THE  PRODUCT  OF  MOLE 

*  FRACTION  AND  CP  OF  THE  SPECIES  DIVIDED  BY  THE  MOLECULAR  WEIGHT  OF 

*  MIXTURE. 

* 

CPS-(A(Q.0)+T(L)*(A(Q,1)+T(L)*(A(Q,2)+T(L)*(A(Q.3)+ 

+  T(L)*A(Q,4)))))*4.186E+07 

CP(L)-CP(L)+CPS*MF(Q,L) 

CPS-CPS/MW(Q) 

H(Q,L)-((B(Q.0)+T(L)*(B(Q,1)+T(L)*(B(Q,2)+T(L)*(B(Q,3)+ 

+  T(L)*B(Q,4) ) ) ) )+HFORM(Q) )/MW(Q)*4 . 186E+07*1000 . 0 

* 

*  COMPUTE  THE  THERMAL  CONDUCTIVITY  OF  THE  MIXTURE 

*  AND  THERMAL  CONDUCTIVITIES  OF  SPECIES  BY  UTILIZATION  OF  THE 

*  EUKEN  FORMULAS.  SEE  BS&L;  EQNS  8.3-15,  8.3-17,  8.3-18,  PP253-260. 

* 

KTS-(CPS+1 . 25*RBAR/MW(Q) )*MUS (Q) 

KT ( L) -KT (L) +MF (Q , L) *KTS/XPHI 

* 

*  COMPUTE  THE  DIFFUSIVITY  OF  SPECIES  BY  CALCULATING 

*  RELATIVE  DIFFERENCES  OF  MOLE  FRACTIONS,  CALCULATING  THE  BINARY 

*  DIFFUSION  COEFFICIENTS,  AND  THEN  EMPLOYING  THESE  IN  A  FORMULA 

*  DERIVED  FROM  THE  STEFAN -MAXWELL  EQUATIONS.  SEE  BS&L; 

*  EQNS  16.4-13,  16.4-15,  16.4-16,  PP508-513,  EQN  18.4-25,  PP569-572. 

* 

IF(DIFSUM . EQ . 0 . 0) THEN 
DIF(Q,L)-0.0 
ELSE 

DIF(Q,L)-(1.0-MF(Q,L))/DIFSUM 
END  IF 
15  CONTINUE 

CP(L)-CP(L)/MWM(L) 

20  CONTINUE 
RETURN 
END 

* 


************************************************************************* 


* 

* 


* 

* 


*  THIS  SUBROUTINE  COMPUTES  THE  MASS  CREATION  TERM  FROM  CHEMICAL  * 

*  REACTION  RATE  EQUATIONS  PROVIDED  BY  THE  USER.  NOTE  THAT  THIS  * 

*  PROCEDURE  IS  CURRENTLY  SET  UP  TO  HANDLE  THE  FOLLOWING  SET  OF  * 

*  REACTIONS  FOR  CHEMICAL  OXYGEN-IODINE  LASERS,  COIL:  * 

*  SEE  V&K;  210-228.  * 

*  * 

*  * 


*********************************************************************** 


* 


* 

* 


1.  02(1D)  +  02(1D) 

2.  02(1S)  +  H20 


-->  02(1S)  +  02(3S) 
-->  02(3S)  +  H20 
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* 

3. 

I2(X)  +  02(1S) 

I  +  I  +  02(3S) 

* 

4. 

I2(X)  +  02(1D) 

I2(STAR)  +  02(3S) 

* 

5. 

12 (X)  +I(STAR) 

-  -> 

I2(STAR)  +  I 

* 

6. 

12 (STAR)  +  02 (ID) 

I  +  I  +  02(3S) 

* 

7. 

12 (STAR)  +  H20 

--> 

12 (X)  +H20 

* 

8. 

I  +  02(1D) 

I (STAR)  +  02 (3S) 

* 

9. 

I(STAR)  +  02(3S) 

I  +  02(1D) 

* 

10. 

I (STAR)  +  02 (ID) 

--> 

I  +  02(1S) 

* 

* 

* 

11. 

I (STAR)  +  H20 

I  +  H20 

SUBROUTINE  RATE(W,R,T,WD,MW) 

INTEGER  L.N.Q.NS 
PARAMETER  (NS-10,N-32) 

REAL*8  W(NS,0:N) .R(0:N) ,T(0:N) ,WD(NS,0;N) ,MW(NS) , F(NS) 

+  .K(ll) 

REAL*8  NAV 

PARAMETER  (NAV-6 . 02252E+23) 

* 

*  COMPUTE  THE  RATE  COEFFICIENTS. 

* 

K(1)-2.7E-17*NAV 

K(2)-5.5E-12*NAV 

K(3)-4.0E-12*NAV 

K(4)-7.0E-15*NAV 

K(5)-3.6E-11*NAV 

K(6)-3.0E-10*NAV 

K(7)-3.0E-10*NAV 

K(8)-7.6E-11*NAV 

K(9)-7.1E-12*NAV 

K(10)-1.1E-13*NAV 

K(11)-2.0E-12*NAV 

* 

*  START  THE  LOOP  FOR  Y. 

* 

DO  10  L-0,N 

* 

*  COMPUTE  THE  TEMPERATURE  DEPENDENT  RATE  COEFFICIENTS. 

* 

K(1)-9.5E-28*T(L)**3.8*EXP(700.0/T(L))*NAV 

K(5)-1.4E-13*EXP(1660/T(L))*NAV 

K(8)-2.33E-8/T(L)*NAV 

K(9)-3.lE-8/T(L)*EXP(-403/T(L))*NAV 

K(10)-4.0E-24*T(L)**3.8*EXP(700.0/T(L))*NAV 

* 

*  COMPUTE  THE  MOLE-MASS  FRACTIONS  FOR  EACH  SPECIES. 

* 

DO  5  0-1, NS 

F(Q)-W(Q.L)/MW(Q) 

5  CONTINUE 

WD(1,L)-(K(1)*F(2)**2+K(2)*F(3)*F(9)+K(3)*F(4)*F(3)+ 

+  K(4)*F(4)*F(2)+K(6)*F(5)*F(2)+K(8)*F(6)*F(2)- 

+  K(9)*F(7)*F(1))*MW(1)*R(L)**2 

WD(2,L)-(-2.0*K(1)*F(2)**2-K(4)*F(4)*F(2)-K(6)*F(5)*F(2)- 
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+  K(8)*F(6)*F(2)+K(9)*F(7)*F{1)-K(10)*F(7)*F(2)) 

+  *MW(2)*R(L)**2 

WD(3,L)-(K(1)*F(2)**2-K(2)*F(3)*F(9)-K(3)*F(4)*F(3)+ 

+  K(10)*F(7)*F(2))*MW(3)*R(L)**2 

WD(4,L)-(-K(3)*F(4)*F(3)-K(4)*F(4)*F(2)-K(5)*F(4)*F(7)+ 
+  K(7)*F(5)*F(9))*MW(4)*R(L)**2 

WD(5,L)-(K(4)*F(4)*F(2)+K(5)*F(4)*F(7)-K(6)*F(5)*F(2)- 
+  K(7)*F(5)*F(9))*MW(5)*R(L)**2 

WD(6,L)-(2.0*K(3)*F(4)*F(3)+K(5)*F(4)*F(7)+ 

+  2.0*K(6)*F(5)*F(2)-K(8)*F(6)*F(2)+K(9)*F(7)*F(1)+ 

+  K(10)*F(7)*F(2)+K(11)*F(7)*F(9))*MW(6)*R(L)**2 

WD(7,L)-(-K(5)*F(4)*F(7)+K(8)*F(6)*F(2)-K(9)*F(7)*F(l)- 
+  K(10)*F(7)*F(2)-K(11)*F(7)*F(9))*MW(7)*R(L)**2 

WD(8.L)-0.0 
WD(9,L)-0.0 
WD(10,L)-0.0 
10  CONTINUE 
RETURN 
END 
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